Page 1
Gianluigi Rozza, Martin Hess, Giovanni Stabile, Marco Tezzele, and Francesco Ballarin Basic Ideas and Tools for Projection-Based Model Reduction of Parametric Partial Differential Equations Abstract: We provide first the functional analysis background required for re- duced order modeling and present the underlying concepts of reduced basis model reduction. The projection-based model reduction framework under affinity assumptions, offline-online decomposition and error estimation is introduced. Several tools for geometry parametrizations, such as free form deformation, ra- dial basis function interpolation and inverse distance weighting interpolation are explained. The empirical interpolation method is introduced as a general tool to deal with non-affine parameter dependency and non-linear problems. The discrete and matrix versions of the empirical interpolation are considered as well. Active subspaces properties are discussed to reduce high-dimensional parameter spaces as a pre-processing step. Several examples illustrate the methodologies. Keywords: Reduced Basis Method, Error Estimation, Free Form Deformation, Radial Basis Function Interpolation, Inverse Distance Weighting Interpolation, Empirical Interpolation Method, Active Subspaces Introduction Parametric model order reduction techniques have been developed in recent decades to deal with increasingly complex computational tasks. The ability to compute how quantities of interest change with respect to parameter variations Gianluigi Rozza, Mathematics Area, SISSA mathLab, Trieste, Italy, gianluigi.rozza@sissa.it Martin Hess, Mathematics Area, SISSA mathLab, Trieste, Italy, martin.hess@sissa.it Giovanni Stabile, Mathematics Area, SISSA mathLab, Trieste, Italy, giovanni.stabile@sissa.it Marco Tezzele, Mathematics Area, SISSA mathLab, Trieste, Italy, marco.tezzele@sissa.it Francesco Ballarin, Mathematics Area, SISSA mathLab, Trieste, Italy, francesco.ballarin@sissa.it arXiv:1911.08954v1 [math.NA] 20 Nov 2019
Page 2
2 G. Rozza et al. provides insight and understanding, which is vital in all areas of science and engineering. Model reduction thus allows to deal with optimization or inverse problems of a whole new scale. Each chapter of the handbook gives an in-depth view of a model order reduction (MOR, or reduced order modeling (ROM)) method, a particular application area, analytical, numerical or technical aspects of software frameworks for model reduction. There exists a large number of MOR techniques used in many areas of sci- ence and engineering to improve computational performances and contain costs in a repetitive computational environment such as many-query and real-time computing [Schilders et al.(2008)Schilders, van der Vorst, and Rommes]. We as- sume a given parametrized partial differential equation (PDE) as starting point of the model reduction procedure. Typical parameters of interest are material coefficients, corresponding to physical qualities of the media which constitute the domain where the PDE is solved. Also a variable geometry can be of special interest in a task to find the optimal device configuration. Physical states such as the temperature might be considered an input parameter. It is a task of the mathematical modeling to identify the parameters of interest and how they enter the PDE. Once a parametrized model is identified, the MOR techniques described in this and the following chapters can be used either in a ‘black-box’ fashion (non-intrusive way) or by intrusive means, which will be explained in detail, whenever this is necessary. The particular numerical method to solve a PDE is most often not relevant to the model reduction procedure. We will therefore assume there is a numerical method available, which solves the problem to any required accuracy, and move seamlessly from the continuous form to the discretized form. This chapter covers briefly the functional analysis framework relevant to many, but not all, MOR methods. Presented is the starting point of PDE-oriented MOR techniques, which can be found in the various chapters of the handbook. In particular, section 1 provides what is needed for the projection-based ROM. Starting from the setting of the classical Lax-Milgram theorem for elliptic PDEs in subsection 1.1 and subsection 1.2, a numerical discretisation is introduced in subsection 1.2.1. Due to brevity of representation, many concepts of functional analysis and theory of PDEs are only ouched upon. Many references to the literature for further reading are given. Projection-based ROM is presented in subsection 1.3, with the following topics covered in detail: proper orthogonal decomposition in subsection 1.3.1, the greedy algorithm in subsection 1.3.2, the projection framework in subsection 1.3.3, affine parameter dependency in subsection 1.3.4, the offline-online decomposition in subsection 1.3.6 and basic error estimation in subsection 1.4.
Page 3
Basic Ideas and Tools for MOR of Parametric PDEs 3 Section 2 introduces efficient techniques for geometric parametrizations, arising from a reference domain approach, such as free form deformation in subsection 2.1, radial basis function interpolation in subsection 2.2 and inverse distance weighting in subsection 2.3. A widely used method to generate an approximate affine parameter depen- dency is the empirical interpolation method. The original empirical interpolation method is presented in section 3 as well as the discrete empirical interpolation method in subsection 3.3 and further options in subsection 3.4. Several numerical examples show the use of the empirical interpolation method in subsection 3.5. Section 4 introduces active subspaces as a pre-processing step to reduce the parameter space dimension. Corresponding examples are provided in sub- section 4.3 and also nonlinear dimensionality reduction is briefly discussed in subsection 4.5. A brief conclusion and outlook at the handbook is given in section 5. 1 Basic Notions and Tools We briefly cover a few main results of linear functional analysis and the analysis of partial differential equations (PDEs). This material serves as a reminder of the underlying concepts of model reduction but can not replace a textbook on these subjects. For a more thorough background, we refer to the literature on functional analysis [Ciarlet(2014), Yosida(1995)], partial differential equa- tions [Adams(1975), Evans(1998), Renardy and Rogers(2004), Rudin(1976)], and numerical methods [Ainsworth and Oden(1997), Babuska(1970/71), Ciarlet(2002), Graetsch and Bathe(2005), Quarteroni(2017), Verfuerth(2013)]. 1.1 Parametrized Partial Differential Equations Let Ω ⊂R𝑑denote a spatial domain in 𝑑= 1, 2 or 3 dimensions with boundary 𝜕Ω. A Dirichlet boundary Γ𝐷⊂𝜕Ω is given, where essential boundary conditions on the field of interest are prescribed. Introduce a Hilbert space 𝑉(Ω) equipped with inner product (·, ·)𝑉and induced norm ‖ · ‖𝑉. A Hilbert space 𝑉(Ω) is a function space, i.e., a function 𝑢∈𝑉(Ω) is seen as a point in the vector space 𝑉, as is common in functional analysis. Each 𝑢∈𝑉(Ω) defines a mapping 𝑥∈Ω ↦→𝑢(𝑥) ∈R or 𝑥∈Ω ↦→𝑢(𝑥) ∈C, depending on whether a real or complex Hilbert space is considered. In many applications, 𝑉is a subset of the Sobolev space 𝐻1(Ω) as 𝑉(Ω) = {𝑣∈𝐻1(Ω): 𝑣|Γ𝐷= 0}. Vector-valued
Page 4
4 G. Rozza et al. Hilbert spaces can be constructed using the Cartesian product of 𝑉(Ω). Given a parameter domain 𝒫⊂R𝑝, a particular parameter point is denoted by the 𝑝-tuple 𝜇= (𝜇1, 𝜇2, … , 𝜇𝑝). The set of all linear and continous forms on 𝑉defines the dual space 𝑉′ and let 𝐿∈ℒ(𝑉, 𝑉′) denote a linear differential operator. A field variable 𝑢∈𝑉: Ω →R is defined implicitly as the solution to a parametrized linear partial differential equation (PDE) through the operator 𝐿: 𝑉× 𝒫→𝑉′ with 𝐿(·; 𝜇) ∈ℒ(𝑉, 𝑉′) and load vector 𝑓𝐿(𝜇) ∈𝑉′ for each fixed 𝜇, as 𝐿(𝑢; 𝜇) = 𝑓𝐿(𝜇). (1) As in the case of function spaces, operators between function spaces form vector spaces themselves, such as 𝐿(·; 𝜇) ∈ℒ(𝑉, 𝑉′), with ℒ(𝑉, 𝑉′) being the space of operators mapping from the vector space 𝑉to 𝑉′. Typical examples of scalar-valued linear PDEs are the Poisson equation, Heat equation or Wave equation, while typical examples of vector-valued linear PDEs are Maxwells equations or Stokes equations. The nonlinear case will be addressed in various chapters as well: examples of nonlinear PDEs are the Navier-Stokes system or the equations describing nonlinear elasticity. 1.2 Parametrized Variational Formulation The variational or weak form of a parametrized linear PDE in the continuous setting is given as 𝑎(𝑢(𝜇), 𝑣; 𝜇) = 𝑓(𝑣; 𝜇) ∀𝑣∈𝑉, (2) with bilinear form 𝑎: 𝑉× 𝑉× 𝒫→R and linear form 𝑓: 𝑉× 𝒫→R. In many application scenarios, a particular output of interest is sought, given by the linear form 𝑙: 𝑉× 𝒫→R as 𝑠(𝜇) = 𝑙(𝑢(𝜇); 𝜇). (3) In the case that 𝑎(·, ·; 𝜇) is symmetric and 𝑙= 𝑓, the problem is called compliant. For each 𝜇∈𝒫assume coercivity and continuity of the bilinear form 𝑎(·, ·; 𝜇), i.e., 𝑎(𝑤, 𝑤; 𝜇) ≥ 𝛼(𝜇)‖𝑤‖2 𝑉, (4) 𝑎(𝑤, 𝑣; 𝜇) ≤ 𝛾(𝜇)‖𝑤‖𝑉‖𝑣‖𝑉, (5)
Page 5
Basic Ideas and Tools for MOR of Parametric PDEs 5 and continuity of the linear form 𝑓(·; 𝜇), 𝑓(𝑤; 𝜇) ≤𝛿(𝜇)‖𝑤‖𝑉, (6) with parameter-independent bounds, which satisfy 0 < 𝛼≤𝛼(𝜇), 𝛾(𝜇) ≤𝛾< ∞ and 𝛿(𝜇) ≤𝛿< ∞. To do actual computations, the biliner form is discretized into a linear equation. The coercivity property means that the matrix discretizing the bilinear form will be positive definite. For fixed parameter the well-posedness of (2) is then established by the Lax-Milgram theorem: Theorem 1.1. Lax-Milgram Theorem Let 𝑎: 𝑉× 𝑉→R be a continuous and coercive bilinear form over a Hilbert space 𝑉and 𝑓∈𝑉′ a continuous linear form. Then the variational problem 𝑎(𝑢, 𝑣) = 𝑓(𝑣) ∀𝑣∈𝑉, (7) has a unique solution 𝑢∈𝑉and it holds ‖𝑢‖𝑉≤1 𝛼‖𝑓‖𝑉, (8) with the coercivity constant 𝛼> 0 of the bilinear form. Thus, in the parametric setting, the 𝜇-dependence also carries over to the coercivity constant as 𝛼= 𝛼(𝜇). The function space in which the field variable resides is called the ansatz space, while the second function space is called the test space, i.e., where a test function 𝑣resides. If the test space is distinct from the ansatz space then the bilinear form is defined over 𝑎: 𝑉× 𝑊× 𝒫→R for 𝑉and 𝑊Hilbert spaces. With 𝑓∈𝑊′ and for fixed 𝜇, the well-posedness is then established through the Banach-Nečas-Babuška theorem. Theorem 1.2. Banach-Nečas-Babuška Theorem Let 𝑉and 𝑊denote Hilbert spaces, 𝑎: 𝑉× 𝑊→R a continuous bilinear form and 𝑓∈𝑊′. Then the variational problem 𝑎(𝑢, 𝑣) = 𝑓(𝑣) ∀𝑣∈𝑊, (9) has a unique solution if and only if (i) the inf-sup condition holds, i.e, ∃𝛽> 0 , s.t., 𝛽≤ inf 𝑣∈𝑉∖{0} sup 𝑤∈𝑊∖{0} 𝑎(𝑣, 𝑤) ‖𝑣‖𝑉‖𝑤‖𝑊 ,
Page 6
6 G. Rozza et al. (ii) ∀𝑤∈𝑊: {𝑎(𝑣, 𝑤) = 0 ∀𝑣∈𝑉} =⇒𝑤= 0. 1.2.1 Discretized Parametrized Variational Formulation The method of weighted residuals is used to cast (1) into a discrete variational formulation. Given the linear PDE 𝐿(𝑢; 𝜇) = 𝑓𝐿(𝜇), consider a discrete, i.e., finite dimensional approximation 𝑢ℎ∈𝑉ℎ⊂𝑉to 𝑢as 𝑢ℎ(𝜇) = 𝑁ℎ ∑︁ 𝑖=1 𝑢(𝑖) ℎ𝜙𝑖. (10) The dimension of 𝑉ℎis 𝑁ℎand the set of ansatz functions 𝜙𝑖(x) : Ω → R belong to 𝑉. The 𝑢(𝑖) ℎ are scalar coefficients such that the vector uℎ= (𝑢(1) ℎ, … , 𝑢(𝑁ℎ) ℎ )𝑇∈R𝑁ℎis the coordinate representation of 𝑢ℎin the basis {𝜙𝑖} of 𝑉ℎ. A conforming discretizations is considered, i.e., 𝑉ℎ⊂𝑉holds. Plugging (10) into (1) yields the discrete residual 𝑅(𝑢ℎ(𝜇)) = 𝐿(𝑢ℎ(𝜇); 𝜇) − 𝑓𝐿(𝜇) ∈𝑉′. To compute the scalar coefficients 𝑢(𝑖) ℎ, Galerkin orthogonality is invoked, as 0 = (𝜙𝑗, 𝑅)(𝑉,𝑉′), 𝑗= 1 … 𝑁ℎ, (11) where (·, ·)(𝑉,𝑉′) is the duality pairing between 𝑉and 𝑉′. In short, Galerkin orthogonality means that the test space is orthogonal to the residual. In Ritz-Galerkin methods, the residual is tested against the same set of functions as the ansatz functions. If test space and trial space are different, one speaks of a Petrov-Galerkin method. Numerous discretization methods can be understood in terms of the method of weighted residuals. They are distinguished by the particular choice of trial and test space. The well-posedness of the discrete setting follows the presentation of the continuous setting, by casting the equations and properties over 𝑉ℎinstead of 𝑉. The weak form in the discrete setting is given as 𝑎(𝑢ℎ(𝜇), 𝑣ℎ; 𝜇) = 𝑓(𝑣ℎ; 𝜇) ∀𝑣ℎ∈𝑉ℎ, (12) with bilinear form 𝑎: 𝑉ℎ× 𝑉ℎ× 𝒫→R and linear form 𝑓: 𝑉ℎ× 𝒫→R. The discrete bilinear form is then derived from (11) through the integration-by-parts formula and Green’s theorem.
Page 7
Basic Ideas and Tools for MOR of Parametric PDEs 7 Correspondingly, the discrete coercivity constant 𝛼ℎ(𝜇) and the discrete continuity constant 𝛾ℎ(𝜇) are defined 𝛼ℎ(𝜇)
min 𝑤ℎ∈𝑉ℎ 𝑎(𝑤ℎ, 𝑤ℎ; 𝜇) ‖𝑤ℎ‖2 𝑉ℎ , (13) 𝛾ℎ(𝜇)
max 𝑤ℎ∈𝑉ℎmax 𝑣ℎ∈𝑉ℎ 𝑎(𝑤ℎ, 𝑣ℎ; 𝜇) ‖𝑤ℎ‖𝑉ℎ‖𝑣ℎ‖𝑉ℎ . (14) The well-posedness of (2) is then analogously established by the Lax-Milgram theorem and Banach-Nečas-Babuška theorem. Cea’s Lemma is a fundamental result about the approximation quality that can be achieved: Lemma 1.3. Cea’s Lemma Let 𝑎: 𝑉× 𝑉→R be a continuous and coercive bilinear form over a Hilbert space 𝑉and 𝑓∈𝑉′ a continuous linear form. Given a conforming finite-dimensional subspace 𝑉ℎ⊂𝑉, the continuity constant 𝛾and coercivity constant 𝛼of 𝑎(·, ·) it holds for the solution 𝑢ℎto 𝑎(𝑢ℎ, 𝑣ℎ) = 𝑓(𝑣ℎ) ∀𝑣ℎ∈𝑉ℎ, (15) that ‖𝑢−𝑢ℎ‖𝑉≤𝛾 𝛼 inf 𝑣ℎ∈𝑉ℎ‖𝑢−𝑣ℎ‖𝑉. (16) The stiffness matrix Aℎ∈R𝑁ℎ×𝑁ℎassembles the bilinear form entrywise as (Aℎ)𝑖𝑗= 𝑎(𝜙𝑗, 𝜙𝑖). The load vector fℎ∈R𝑁ℎis assembled entrywise as (fℎ)𝑖= 𝑓(𝜙𝑖) and the solution vector is denoted uℎwith coefficients 𝑢(𝑗) ℎ. Then solving (12) amounts to solving the linear system Aℎuℎ= fℎ. (17) The most common discretization method is the finite element method (FEM) [Boffi et al.(2013)Boffi, Brezzi, and Fortin], besides finite difference [Smith(1985)], discontinuous Galerkin [Arnold et al.(2002)Arnold, Brezzi, Cockburn, and Marini], finite volume [Eymard et al.(2000)Eymard, Gallouët, and Herbin] and spectral element method [Canuto et al.(2006)Canuto, Hussaini, Quarteroni, and Zhang]. 1.3 Model Reduction Basic Concepts A wide variety of reduced order modeling methods exist today, thanks to large research efforts in the last decades. Reduced basis model order reduction is a
Page 8
8 G. Rozza et al. projection-based MOR method and also shares many features with other MOR methods, so that the topics mentioned here will occur throughout the handbook. Two common algorithms for the generation of a projection space, the proper orthogonal decomposition and the greedy algorithm, are presented first. 1.3.1 Proper Orthogonal Decomposition Assume a sampled set of high-fidelity solutions {𝑢ℎ(𝜇𝑖), 𝑖= 1, … , 𝑁𝑚𝑎𝑥}, i.e., solutions to (12) or (17), respectively. The discrete solution vectors are stored column-wise in a snapshot matrix S ∈R𝑁ℎ×𝑁𝑚𝑎𝑥. The proper orthogonal de- composition (POD) compresses the data stored in S by computing an orthogonal matrix V, which is a best approximation in the least square sense to S. In particular, the POD solution of size 𝑁is the solution to min V∈R𝑁ℎ×𝑁‖S −VV𝑇S‖𝐹, (18) subject to V𝑇V = I𝑁×𝑁, (19) with ‖ · ‖𝐹the Frobenius norm and I𝑁×𝑁the identity matrix. There exists a solution to (18) - (19) according to the Eckardt-Young-Mirsky theorem [Eckart and Young(1936)], which can be computed with the singular value decomposition (SVD) as S = UΣZ, (20) with orthogonal matrix U ∈R𝑁ℎ×𝑁ℎ, rectangular diagonal matrix Σ ∈ R𝑁ℎ×𝑁𝑚𝑎𝑥and orthogonal matrix Z ∈R𝑁𝑚𝑎𝑥×𝑁𝑚𝑎𝑥. The solution V is com- posed of the first 𝑁column vectors of U. They are also called the POD modes. The diagonal entries {𝜎𝑖, 𝑖= 1, … , min(𝑁ℎ, 𝑁𝑚𝑎𝑥)} of Σ are non-negative and called singular values. It holds that min V∈R𝑁ℎ×𝑁‖S −VV𝑇S‖𝐹= min(𝑁ℎ,𝑁𝑚𝑎𝑥) ∑︁ 𝑖=𝑁+1 𝜎𝑖. (21) Thus, the neglected singular values give an indication of the approximate truncation error. In practise, a high tolerance threshold like 99% or 99.99% is chosen and 𝑁is determined so that the sum of the first 𝑁singular values reaches this percentage of the sum of all singular values. In many applications, an exponential singular value decay can be observed, which allows to reach the tolerance with a few POD modes.
Page 9
Basic Ideas and Tools for MOR of Parametric PDEs 9 1.3.2 Greedy Algorithm The greedy algorithm also computes an orthogonal matrix V ∈R𝑁ℎ×𝑁to serve as a projection operator, just as in the POD case. The greedy algorithm is an iterative procedure, which enriches the snapshot space according to where an error indicator or error estimator Δ attains its maximum. Starting from a field solution at a given initial parameter value, the parameter location is sought, whose field solution is worst approximated with the initial solution. This solution is then computed and appended to the projection matrix to obtain a 2-dimensional projection space. The greedy typically searches for new snapshot solutions within a discrete surrogate 𝑃of the parameter space 𝒫. The process is repeated until a given tolerance on the error estimator is fulfilled. The error estimator is residual-based and estimates the error between a reduced order solve for a projection space V and the high-fidelity solution, see subsection 1.4. The greedy algorithm is stated in pseudocode in algorithm 1. Algorithm 1 The greedy algorithm Input: discrete surrogate 𝑃of parameter space 𝒫, approximation tolerance tol, initial parameter 𝜇1 Output: projection matrix V 𝑁= 1 V1 = uℎ(𝜇1) ‖uℎ(𝜇1)‖ while max𝜇∈𝑃Δ(𝜇) > tol do 𝑁= 𝑁+ 1 𝜇𝑁= arg max 𝜇∈𝑃 Δ(𝜇) solve (17) at 𝜇𝑁for uℎ(𝜇𝑁) orthonormalize uℎ(𝜇𝑁) with respect to V𝑁−1 to obtain 𝜁𝑁 append 𝜁𝑁to V𝑁−1 to obtain V𝑁 end while set V = V𝑁 1.3.3 Reduced Order System Starting from the discrete high-fidelity formulation (12), another Galerkin projec- tion is invoked to arrive at the reduced order formulation. Assume a projection space 𝑉𝑁is then determined through either a proper orthogonal decomposition
Page 10
10 G. Rozza et al. (POD) or the greedy sampling, with V ∈R𝑁ℎ×𝑁denoting a discrete basis of 𝑉𝑁. Thus 𝑉𝑁⊂𝑉ℎand dim 𝑉𝑁= 𝑁. The reduced order variational formulation is to determine 𝑢𝑁(𝜇) ∈𝑉𝑁, such that 𝑎(𝑢𝑁(𝜇), 𝑣𝑁; 𝜇) = 𝑓(𝑣𝑁; 𝜇) ∀𝑣𝑁∈𝑉𝑁. (22) Eq. (17) is then projected onto the reduced order space as V𝑇AℎVu𝑁= V𝑇fℎ. (23) The reduced system matrix A𝑁= V𝑇AℎV is then a dense matrix of small size 𝑁× 𝑁as depicted in (24): [︀ A𝑁 ]︀
⎡ ⎣V ⎤ ⎦ 𝑇⎡ ⎣ 𝑎(𝜙1, 𝜙1) … … … … … … … 𝑎(𝜙𝑁ℎ, 𝜙𝑁ℎ) ⎤ ⎦ ⎡ ⎣V ⎤ ⎦. (24) The high-order solution is then approximated as uℎ≈Vu𝑁. (25) 1.3.4 Affine Parameter Dependency Many MOR algorithms rely on an affine parameter dependency, because the affine parameter dependency provides the computational efficiency of the model reduc- tion. Thus, it is a significant advancement from the 2000’s [Rozza et al.(2008)Rozza, Huynh, an over the first use of reduced order models [Almroth et al.(1978)Almroth, Stern, and Brogan, Noor(1982)]. An affine parameter dependency means that the bilinear form can be ex- panded as 𝑎(·, ·; 𝜇) = 𝑄𝑎 ∑︁ 𝑖=1 Θ𝑖 𝑎(𝜇)𝑎𝑖(·, ·), (26) and affine expansions hold as 𝑓(·; 𝜇) = 𝑄𝑓 ∑︁ 𝑖=1 Θ𝑖 𝑓(𝜇)𝑓𝑖(·), (27) 𝑙(·; 𝜇) = 𝑄𝑙 ∑︁ 𝑖=1 Θ𝑖 𝑙(𝜇)𝑙𝑖(·), (28)
Page 11
Basic Ideas and Tools for MOR of Parametric PDEs 11 with scalar-valued functions Θ𝑖 𝑎: 𝒫→R, Θ𝑖 𝑓: 𝒫→R and Θ𝑖 𝑙: 𝒫→R. Correspondingly the linear system (17) can be expanded as (︃𝑄𝑎 ∑︁ 𝑖=1 Θ𝑖 𝑎(𝜇)A𝑖 )︃ uℎ= 𝑄𝑓 ∑︁ 𝑖=1 Θ𝑖 𝑓(𝜇)f𝑖, (29) as well as the reduced order form (23) V𝑇 (︃𝑄𝑎 ∑︁ 𝑖=1 Θ𝑖 𝑎(𝜇)A𝑖 )︃ Vu𝑁= V𝑇 𝑄𝑓 ∑︁ 𝑖=1 Θ𝑖 𝑓(𝜇)f𝑖, (30) (︃𝑄𝑎 ∑︁ 𝑖=1 Θ𝑖 𝑎(𝜇)V𝑇A𝑖V )︃ u𝑁= 𝑄𝑓 ∑︁ 𝑖=1 Θ𝑖 𝑓(𝜇)V𝑇f𝑖. (31) MOR relies on an affine parameter dependency, such that all computations de- pending on the high-order model size can be moved into a parameter-independent offline phase, while having a fast input-output evaluation online. If the problem is not affine, an affine representation can be approximated using a technique such as the empirical interpolation method, see section 3. 1.3.5 Affine shape parametrizations: an example Consider heat conduction in a square domain Ω(𝑥, 𝑦) = [0, 1]2. On the left side 𝑥= 0, inhomogeneous Neumann conditions, i.e., a non-zero heat flux is imposed and on the right side 𝑥= 1, homogeneous Dirichlet conditions, i.e., zero temperature is imposed. The top and bottom side, homogeneous Neumann conditions, i.e., a zero heat flux is imposed. Consider two different media with different conductivities 𝜎1 and 𝜎2 occupying the subdomains Ω1(𝜇) = [0, 𝜇]×[0, 1] and Ω2(𝜇) = [𝜇, 1] × [0, 1], for 𝜇∈𝒫= (0, 1), as shown in Figure 1. For the sake of clarity, in the rest of this section we identify the 1-dimensional parameter vector 𝜇with its (only) component 𝜇, thus dropping the bold notation from the symbol. Choosing 𝜇= 0.5 as the reference configuration, there exist affine transfor- mations from the reference domain to the actual domain. It holds 𝑇1 : Ω1(𝜇) →Ω1(𝜇) : (𝑥, 𝑦) ↦→(2𝜇𝑥, 𝑦), (32) 𝑇2 : Ω2(𝜇) →Ω2(𝜇) : (𝑥, 𝑦) ↦→((2 −2𝜇)𝑥, 𝑦) + (2𝜇−1, 0). (33) In general, an affine transformation of a subdomain can be expressed as 𝑇𝑘: Ω𝑘(𝜇) →Ω𝑘(𝜇) : x ↦→𝐺𝑘(𝜇)x + 𝐷𝑘(𝜇), (34)
Page 12
12 G. Rozza et al. Fig. 1: The computational domain is subdivided into two domains Ω = Ω1 ∪Ω2, depending on the parameter 𝜇. Shown here for 𝜇= 0.5. with x ∈R𝑑, 𝐺𝑘∈R𝑑×𝑑and 𝐷𝑘∈R𝑑in 𝑑= 2, 3 spatial dimensions. Thus, the bilinear form 𝑎(𝑢, 𝑣; 𝜇) = ∫︁ Ω1(𝜇) 𝜎1∇𝑢· ∇𝑣𝑑x + ∫︁ Ω2(𝜇) 𝜎2∇𝑢· ∇𝑣𝑑x (35) can be mapped to the reference domain with the inverse affine transformation 𝑇−1 𝑘 : Ω𝑘(𝜇) →Ω𝑘(𝜇) : x ↦→𝐺−1 𝑘(𝜇)x −𝐺−1 𝑘(𝜇)𝐷𝑘(𝜇), (36) and integration by substitution as 𝑎(𝑢, 𝑣; 𝜇) = ∫︁ Ω1(𝜇) 𝜎1(∇𝑢𝐺−1 1 (𝜇)) · (𝐺−𝑇 1 (𝜇)∇𝑣) det(𝐺1(𝜇)) 𝑑x (37) + ∫︁ Ω2(𝜇) 𝜎2(∇𝑢𝐺−1 2 (𝜇)) · (𝐺−𝑇 2 (𝜇)∇𝑣) det(𝐺2(𝜇)) 𝑑x, (38) which establishes the affine parameter dependency (26) by computing Θ𝑖 𝑎(𝜇) from the coefficients of 𝐺1 and 𝐺2 [Rozza et al.(2008)Rozza, Huynh, and Patera, Rozza et al.(2009)Rozza, Huynh, Nguyen, and Patera, Chinesta et al.(2017)Chinesta, Huert It is:
Page 13
Basic Ideas and Tools for MOR of Parametric PDEs 13 ∫︁ Ω1(𝜇) 𝜎1(∇𝑢𝐺−1 1 (𝜇)) · (𝐺−𝑇 1 (𝜇)∇𝑣) det(𝐺1(𝜇)) 𝑑x (39)
∫︁ Ω1(𝜇) 𝜎1((2𝜇)−1𝜕𝑥𝑢, 𝜕𝑦𝑢) · ((2𝜇)−1𝜕𝑥𝑣, 𝜕𝑦𝑣)2𝜇𝑑x (40) = (2𝜇)−1 ∫︁ Ω1(𝜇) 𝜎1(𝜕𝑥𝑢)(𝜕𝑥𝑣) 𝑑x + 2𝜇 ∫︁ Ω1(𝜇) 𝜎1(𝜕𝑦𝑢)(𝜕𝑦𝑣) 𝑑x, (41) and ∫︁ Ω2(𝜇) 𝜎2(∇𝑢𝐺−1 2 (𝜇)) · (𝐺−𝑇 2 (𝜇)∇𝑣) det(𝐺2(𝜇)) 𝑑x (42)
∫︁ Ω2(𝜇) 𝜎2((2 −2𝜇)−1𝜕𝑥𝑢, 𝜕𝑦𝑢) · ((2 −2𝜇)−1𝜕𝑥𝑣, 𝜕𝑦𝑣)(2 −2𝜇) 𝑑x (43) = (2 −2𝜇)−1 ∫︁ Ω2(𝜇) 𝜎2(𝜕𝑥𝑢)(𝜕𝑥𝑣) 𝑑x + (2 −2𝜇) ∫︁ Ω2(𝜇) 𝜎2(𝜕𝑦𝑢)(𝜕𝑦𝑣) 𝑑x, (44) which establishes the affine form (26) with 𝑄𝑎= 4 and Θ1 𝑎(𝜇) = (2𝜇)−1, (45) Θ2 𝑎(𝜇) = 2𝜇, (46) Θ3 𝑎(𝜇) = (2 −2𝜇)−1, (47) Θ4 𝑎(𝜇) = 2 −2𝜇, (48) and
Page 14
14 G. Rozza et al. 𝑎1(·, ·) = ∫︁ Ω1(𝜇) 𝜎1(𝜕𝑥𝑢)(𝜕𝑥𝑣) 𝑑x, (49) 𝑎2(·, ·) = ∫︁ Ω1(𝜇) 𝜎1(𝜕𝑦𝑢)(𝜕𝑦𝑣) 𝑑x, (50) 𝑎3(·, ·) = ∫︁ Ω2(𝜇) 𝜎2(𝜕𝑥𝑢)(𝜕𝑥𝑣) 𝑑x, (51) 𝑎4(·, ·) = ∫︁ Ω2(𝜇) 𝜎2(𝜕𝑦𝑢)(𝜕𝑦𝑣) 𝑑x. (52) The second and fourth term can be further simplified to a term depending on 2𝜇and a 𝜇-independent term, but in this case it still leaves 𝑄𝑎= 4 terms. In some cases the number of affine terms can be automatically reduced further using symbolic computations. 1.3.6 Offline-Online Decomposition The offline-online decomposition enables the computational speed-up of the ROM approach in many-query scenarios. It is also known as the offline-online paradigm, which assumes that a compute-intensive offline phase can be performed on a supercomputer, which generates all quantities depending on the large discretization size 𝑁ℎ. Once completed, a reduced order solve, i.e., an online solve for a new parameter of interest can be performed with computational cost independent of the large discretization size 𝑁ℎ. The online phase can thus be performed even on mobile and embedded devices, see Figure 2. If a supercomputer is not available, this can be relaxed however. There exist heuristic algorithms to make also the offline phase feasible on a common workstation, such that a typical scenario would be that the offline phase runs overnight and a reduced model is available the next morning. Noticing that the terms V𝑇A𝑖V and V𝑇f𝑖in (31) are parameter-independent, they can be precomputed, prior to any ROM parameter sweep. This will store small-sized dense matrices of dimension 𝑁× 𝑁. Once a reduced order solution u𝑁is desired for a given parameter 𝜇, the sum given in (31) is formed and solved for u𝑁. Since this is the same as solving (23), the reduced order approximation is then available as uℎ≈Vu𝑁, see (25).
Page 15
Basic Ideas and Tools for MOR of Parametric PDEs 15 Fig. 2: Offline-online paradigm. The complex high fidelity simulations are carried out in high performance (HPC) clusters for given preselected parameters. The solution snapshots can be stored and the ROM trained. Then in the offline phase the ROM provides approximated solutions at new untried parameters in real time on simple portable devices. 1.4 Error bounds In this section we develop effective and reliable a posteriori error estimators for the field variable or an output of interest. The use of such error bounds drives the construction of the reduced basis during the offline stage, thanks to the so-called greedy algorithm. Moreover, during the online stage, such bounds provide a certified accuracy of the proposed reduced order model. Following [Rozza et al.(2008)Rozza, Huynh, and Patera], we introduce residual-based a posteriori error estimation for the elliptic case. From (12) and (22) it follows that the error 𝑒(𝜇) = 𝑢ℎ(𝜇) −𝑢𝑁(𝜇) satisfies 𝑎(𝑒(𝜇), 𝑣ℎ; 𝜇) = 𝑟(𝑣ℎ; 𝜇) ∀𝑣ℎ∈𝑉ℎ. (53) where the residual 𝑟(·; 𝜇) ∈𝑉′ ℎis defined as 𝑟(𝑣ℎ; 𝜇) = 𝑓(𝑣ℎ; 𝜇) −𝑎(𝑢𝑁(𝜇), 𝑣ℎ; 𝜇) ∀𝑣ℎ∈𝑉ℎ. (54) The following theorem further characterizes the relation between error and residual: Theorem 1.4. Under compliance assumptions, the following inequalities hold ‖𝑒(𝜇)‖𝜇= ‖𝑢ℎ(𝜇) −𝑢𝑁(𝜇)‖𝜇≤Δ𝑒𝑛(𝜇) = ‖𝑟(·; 𝜇)‖𝑉′ ℎ √︀ 𝛼ℎ(𝜇) , (55) 0 ≤𝑠ℎ(𝜇) −𝑠𝑁(𝜇) ≤Δ𝑠(𝜇) = ‖𝑟(·; 𝜇)‖2 𝑉′ ℎ 𝛼ℎ(𝜇) , (56) where ‖𝑣‖2 𝜇= 𝑎(𝑣, 𝑣; 𝜇) defines an equivalent norm to ‖𝑣‖𝑉ℎ. Proof. ‖·‖𝜇defines an equivalent norm thanks to symmetry, continuity and coercivity of 𝑎(·, ·; 𝜇).
Page 16
16 G. Rozza et al. Since 𝑒(𝜇) ∈𝑉ℎ, from (53) with 𝑣ℎ= 𝑒(𝜇) it follows ‖𝑒(𝜇)‖2 𝜇= 𝑎(𝑒(𝜇), 𝑒(𝜇); 𝜇) = 𝑟(𝑒(𝜇); 𝜇) ≤‖𝑟(·; 𝜇)‖𝑉′ ℎ‖𝑒(𝜇)‖𝑉ℎ, the last inequality being due to the definition of the norm in 𝑉′ ℎ. Furthermore, due to coercivity, it holds ‖𝑒(𝜇)‖2 𝜇= 𝑎(𝑒(𝜇), 𝑒(𝜇); 𝜇) ≥𝛼(𝜇) ‖𝑒(𝜇)‖2 𝑉ℎ. Combining these two results yields (55). Furthermore, since 𝑙= 𝑓are linear forms, 𝑠ℎ(𝜇) −𝑠𝑁(𝜇) = 𝑙(𝑒(𝜇); 𝜇) = 𝑓(𝑒(𝜇); 𝜇) = 𝑎(𝑢ℎ(𝜇), 𝑒(𝜇); 𝜇) (57) From (53) with 𝑣ℎ:= 𝑣𝑁∈𝑉𝑁and (22) it follows 𝑎(𝑒(𝜇), 𝑣𝑁; 𝜇) = 𝑟(𝑣𝑁(𝜇); 𝜇) = 0. This holds in particular for 𝑣𝑁= 𝑢𝑁(𝜇). Moreover, due to symmetry, 𝑎(𝑢𝑁(𝜇), 𝑒(𝜇); 𝜇) = 0 as well. Thus, 𝑎(𝑢ℎ(𝜇), 𝑒(𝜇); 𝜇) = 𝑎(𝑒(𝜇), 𝑒(𝜇); 𝜇) in (57), and we conclude that 𝑠ℎ(𝜇) −𝑠𝑁(𝜇) = ‖𝑒(𝜇)‖2 𝜇. (58) The upper bound in (56) is then a consequence of (55), while the lower bound trivially holds as the right-hand side of (58) is a non-negative quantity. Offline-online decomposition is usually solicited for the a posteriori error bounds introduced by the previous theorem, for the sake of a fast computation of the right-hand side of (55)-(56). This requires the efficient evaluation of both the numerator (dual norm of the residual) and the denominator (parametrized coercivity constant). The Riesz representation theorem is employed to define the unique ^𝑟(𝜇) ∈𝑉ℎsuch that (^𝑟(𝜇), 𝑣ℎ)𝑉ℎ= 𝑟(𝑣ℎ; 𝜇), ∀𝑣ℎ∈𝑉ℎ. (59) Under affine separability assumptions (26)-(28), it holds 𝑟(𝑣ℎ; 𝜇) = 𝑄𝑓 ∑︁ 𝑖=1 Θ𝑖 𝑓(𝜇)𝑓𝑖(𝑣ℎ) − 𝑁 ∑︁ 𝑛=1 u𝑁𝑛 𝑄𝑎 ∑︁ 𝑖=1 Θ𝑖 𝑎(𝜇)𝑎𝑖(𝜁𝑛, 𝑣ℎ), ∀𝑣ℎ∈𝑉ℎ;
Page 17
Basic Ideas and Tools for MOR of Parametric PDEs 17 so that an affine expansion with 𝑄𝑓+ 𝑁𝑄𝑎terms is obtained for 𝑟(·; 𝜇). Riesz representation is then invoked for 𝑟1(𝑣ℎ; 𝜇) = 𝑓1(𝑣ℎ), … , 𝑟𝑄𝑓(𝑣ℎ; 𝜇) = 𝑓𝑄𝑓(𝑣ℎ), 𝑟𝑄𝑓+1(𝑣ℎ; 𝜇) = 𝑎1(𝜁1, 𝑣ℎ), … , 𝑟𝑄𝑓+𝑄𝑎(𝑣ℎ; 𝜇) = 𝑎𝑄𝑎(𝜁1, 𝑣ℎ), … 𝑟𝑄𝑓+(𝑁−1)𝑄𝑎+1(𝑣ℎ; 𝜇) = 𝑎1(𝜁𝑁, 𝑣ℎ), … , 𝑟𝑄𝑓+𝑁𝑄𝑎(𝑣ℎ; 𝜇) = 𝑎𝑄𝑎(𝜁𝑁, 𝑣ℎ) during the offline stage, storing the corresponding solutions to (59). For what concerns the evaluation of the denominator of (55)-(56), ex- act evaluation of 𝛼(𝜇) is seldom employed. Instead, an offline-online de- composable lower bound is sought. Early proposals on the topic are available in [Veroy et al.(2002)Veroy, Rovas, and Patera, Prud’Homme et al.(2002)Prud’Homme, Rovas Veroy and Patera(2005), Rozza et al.(2008)Rozza, Huynh, and Patera, Canuto et al.(2009)Can In 2007, the successive constraint method (SCM) was devised in [Huynh et al.(2007)Huynh, Rozz based on successive linear programming approximations, and subsequently ex- tended in [Chen et al.(2008)Chen, Hesthaven, Maday, and Rodríguez, Chen et al.(2009)Chen, Vallaghé et al.(2011)Vallaghé, Fouquembergh, Le Hyaric, and Prud’Homme, Zhang(2011)]. Alternative methodologies based on interpolation techniques have also appeared in recent years in [Hess et al.(2015)Hess, Grundel, and Benner, Manzoni and Negri(2015), Iapichino et al.(2017)Iapichino, Ulbrich, and Volkwein]. A posteriori error estimation can be derived for more general problems as well (including non-coercive linear, nonlinear or time-dependent prob- lems), through application of the Brezzi-Rappaz-Raviart theory. We refer to [Veroy and Patera(2005), Deparis and Rozza(2009), Yano(2014), Manzoni(2014), Rebollo et al.(2017)Rebollo, Ávila, Mármol, Ballarin, and Rozza] for a few rep- resentative cases. To this end, extensions of SCM are discussed in [Chen et al.(2009)Chen, Hestha Hesthaven et al.(2012)Hesthaven, Stamm, and Zhang, Huynh et al.(2010)Huynh, Knezevic, C Chen(2016)]. 2 Geometrical parametrization for shapes and domains In this section we discuss problems characterized by a geometrical parametriza- tion. In particular, a reference domain approach is discussed, relying on a map that deforms the reference domain into the parametrized one. Indeed, while affine shape parametrization (see section 1.3.5 for an example, and [Rozza et al.(2008)Rozza, Huynh, and Patera] for more details) naturally abides
Page 18
18 G. Rozza et al. by the offline-online separability assumption, it often results in very limited deformation of the reference domain, or strong assumptions on the underlying shape. Let Ω ⊂R𝑑, 𝑑= 2, 3, be the reference domain. Let ℳbe a parametric shape morphing function, that is ℳ(𝑥; 𝜇) : R𝑑→R𝑑 (60) which maps the reference domain Ω into the deformed domain Ω(𝜇) as Ω(𝜇) = ℳ(Ω; 𝜇), where 𝜇∈𝒫represents the vector of the geometrical parameters. This map will change accordingly to the chosen shape morphing technique. The case of Section 1.3.5 is representative of an affine map ℳ(·; 𝜇). Instead, in the following we address more general (not necessarily affine) techniques such as the free form deformation (FFD), the radial basis functions (RBF) interpolation, and the inverse distance weighting (IDW) interpolation. From a practical point of view, we recommend the Python package called PyGeM - Python Geometrical Morphing (see [PyGeM(2017)]), which allows an easy integration with the majority of industrial CAD files and the most common mesh files. 2.1 Free form deformation The free form deformation is a widely used parametrization and morphing technique both in academia and in industry. For the original formulation see [Sederberg and Parry(1986)]. More recent works use FFD coupled with reduced basis methods for shape optimization and de- sign of systems modeled by elliptic PDEs (see [Lassila and Rozza(2010)], [Rozza et al.(2013)Ro and [Sieger et al.(2015)Sieger, Menzel, and Botsch]), in naval engineering for the optimization of the bulbous bow shape of cruise ships (in [Demo et al.(2018a)Demo, Tezzele, G in the context of sailing boats in [Lombardi et al.(2012)Lombardi, Parolini, Quarteroni, and R and in automotive engineering in [Salmoiraghi et al.(2018)Salmoiraghi, Scardigli, Telib, and R FFD can be used both for global and local deformations and it is completely independent to the geometry to morph. It acts through the displacement of a lattice of points, called FFD control points, constructed around the domain of interest. In particular it consists in three different steps as depicted in Figure 3. First the physical domain Ω is mapped to ^Ω, the reference one, through the affine map 𝜓. Then the lattice of control points is constructed and the displacements of these points by the map ^𝑇, is what we call geometrical parameters 𝜇. The deformation is propagated to the entire embedded body usually by using Bern-
Page 19
Basic Ideas and Tools for MOR of Parametric PDEs 19 stein polynomials. Finally through the inverse map 𝜓−1 we return back to the parametric physical space Ω(𝜇). Fig. 3: Scheme of the three maps composing the FFD map ℳ. In particular 𝜓maps the physical space to the reference one, then ^𝑇deforms the entire geometry according to the displacements of the lattice control points, and finally 𝜓−1 maps back the reference domain to the physical one. So, recalling Eq. (60), we have the explicit map ℳfor the FFD, that is the composition of the three maps presented, i.e. ℳ(𝑥, 𝜇) = (𝜓−1 ∘̂︀𝑇∘𝜓)(𝑥, 𝜇) = (61) = 𝜓−1 (︃𝐿 ∑︁ 𝑙=0 𝑀 ∑︁ 𝑚=0 𝑁 ∑︁ 𝑛=0 𝑏𝑙𝑚𝑛(𝜓(𝑥))𝑃0 𝑙𝑚𝑛(𝜇𝑙𝑚𝑛) )︃ ∀𝑥∈Ω, (62) where 𝑏𝑙𝑚𝑛are Bernstein polynomials of degree 𝑙, 𝑚, 𝑛in each direction, respec- tively, 𝑃0 𝑙𝑚𝑛(𝜇𝑙𝑚𝑛) = 𝑃𝑙𝑚𝑛+ 𝜇𝑙𝑚𝑛, with 𝑃𝑙𝑚𝑛representing the coordinates of the control point identified by the three indices 𝑙, 𝑚, 𝑛in the lattice of FFD con- trol points. In an offline-online fashion, for a given 𝑥, terms {𝑏𝑙𝑚𝑛(𝜓(𝑥))}𝑙,𝑚,𝑛 can be precomputed during the offline stage, resulting in an inexpensive linear combination of 𝑥-dependent precomputed quantities and 𝜇-dependent control points locations {𝑃0 𝑙𝑚𝑛(𝜇𝑙𝑚𝑛)}𝑙,𝑚,𝑛. The application of 𝜓−1 does not hinder such offline-online approach as 𝜓is affine. We can notice that the deformation does not depend on the topology of the object to be morphed, so this technique is very versatile and non- intrusive, especially for complex geometries or in industrial contexts (see e.g.
Page 20
20 G. Rozza et al. [Salmoiraghi et al.(2016)Salmoiraghi, Ballarin, Corsi, Mola, Tezzele, and Rozza, Rozza et al.(2018)Rozza, Malik, Demo, Tezzele, Girfoglio, Stabile, and Mola]). In the case where the deformation has to satisfy some constraints, like for example continuity constraints, it is possible to increase the number of control points. Often it is the case where at the interface between the undeformed portion of the geometry and the morphed area the continuity has to be prescribed for physical reasons. As an example, in Figure 4 we present a free form deformation of a bulbous bow, where an STL file of a complete hull is morphed continuously by the displacement of only some control points. Fig. 4: Bulbous bow deformation using FFD. In green the FFD control points defining the morphing. 2.2 Radial basis functions interpolation Radial Basis Functions (RBF) represent a powerful tool for nonlinear multivariate approximation, interpolation between nonconforming meshes ([Deparis et al.(2014)Deparis, Fo and for shape parametrization due to its approximation properties (see [Buhmann(2003)]). A radial basis function is any smooth real-valued function ̃︀𝜙: R𝑑→R such that it exists 𝜙: R+ →R and ̃︀𝜙(𝑥) = 𝜙(‖𝑥‖), where ‖·‖ indicates the Euclidean norm in R𝑑. The most widespread radial basis functions are the following:
Page 21
Basic Ideas and Tools for MOR of Parametric PDEs 21 – Gaussian splines ([Buhmann(2003)]) defined as 𝜙(‖𝑥‖) = 𝑒−‖𝑥‖2/𝑅; – thin plate splines ([Duchon(1977)]) defined as 𝜙(‖𝑥‖) = (︂‖𝑥‖ 𝑅 )︂2 ln (︂‖𝑥‖ 𝑅 )︂ ; – Beckert and Wendland 𝐶2 basis ([Beckert and Wendland(2001)]) defined as 𝜙(‖𝑥‖) = (︂ 1 −‖𝑥‖ 𝑅 )︂4 + (︂ 4‖𝑥‖ 𝑅
- 1 )︂ ; – multi-quadratic biharmonic splines ([Sandwell(1987)]) defined as 𝜙(‖𝑥‖) = √︀ ‖𝑥‖2 + 𝑅2; – inverted multi-quadratic biharmonic splines ([Buhmann(2003)]) defined as 𝜙(‖𝑥‖) = 1 √︀ ‖𝑥‖2 + 𝑅2 ; where 𝑅> 0 is a given radius, and the subscript + indicates the positive part. Following [Morris et al.(2008)Morris, Allen, and Rendall, Manzoni et al.(2012)Manzoni, Q given 𝒩𝐶control points situated on the surface of the body to morph, we can generate a deformation by moving some of these points and imposing the new surface which interpolates them. The displacements of the control points represent the geometrical parameters 𝜇. We can now define the map ℳin Eq. (60) for the RBF interpolation technique, that is ℳ(𝑥; 𝜇) = 𝑞(𝑥; 𝜇) + 𝒩𝐶 ∑︁ 𝑖=1 𝛾𝑖(𝜇) 𝜙(‖𝑥−𝑥𝐶𝑖‖), (63) where 𝑞(𝑥; 𝜇) is a polynomial term, generally of degree 1, 𝛾𝑖(𝜇) is the weight associated to the basis function 𝜙𝑖, and {𝑥𝐶𝑖}𝒩𝐶 𝑖=1 are control points selected by the user (denoted by spherical green markers in Figure 5), and 𝑥∈Ω. We underline that in the three dimensional case (63) has 𝑑× 𝒩𝐶+ 𝑑+ 𝑑2 unknowns, which are 𝑑× 𝒩𝐶for the 𝛾𝑖and 𝑑+ 𝑑2 for the polynomial term 𝑞(𝑥; 𝜇) = 𝑐(𝜇) + Q(𝜇)𝑥. To this end we impose the interpolatory constraint ℳ(𝑥𝐶𝑖; 𝜇) = 𝑦𝐶𝑖(𝜇) ∀𝑖∈{1, … , 𝒩𝐶}, (64)
Page 22
22 G. Rozza et al. Fig. 5: Two different views of the same deformed carotid artery model using the radial basis functions (RBF) interpolation technique. The green dots indicate the RBF control points that define the morphing. The black small points highlight the original undeformed geometry. The occlusion of the two branches is achieved through a displacement along the normal direction with respect to the carotid surface of the control points after the bifurcation. where 𝑦𝐶𝑖are the deformed control points obtained applying the displacement 𝜇to 𝑥𝐶𝑖, in particular 𝑥𝐶= [𝑥𝐶1, … , 𝑥𝐶𝒩𝐶] ∈R𝒩𝐶×𝑑, (65) 𝑦𝐶(𝜇) = [𝑦𝐶1(𝜇), … , 𝑦𝐶𝒩𝐶(𝜇)] ∈R𝒩𝐶×𝑑. (66) For the remaining 𝑑+ 𝑑2 unknowns, due to the presence of the polynomial term, we complete the system with additional constraints that represent the conservation of the total force and momentum (see [Buhmann(2003), Morris et al.(2008)Morris, Allen, and Rendall]), as follows 𝒩𝐶 ∑︁ 𝑖=1 𝛾𝑖(𝜇) = 0, (67) 𝒩𝐶 ∑︁ 𝑖=1 𝛾𝑖(𝜇)[𝑥𝐶𝑖]1 = 0, … 𝒩𝐶 ∑︁ 𝑖=1 𝛾𝑖(𝜇)[𝑥𝐶𝑖]𝑑= 0, (68) where the notation [𝑥]𝑑denotes the 𝑑-th component of the vector 𝑥.
Page 23
Basic Ideas and Tools for MOR of Parametric PDEs 23 Following an offline-online strategy, for a given 𝑥, evaluation of 𝜙(‖𝑥−𝑥𝐶𝑖‖), 𝑖= 1, … , 𝒩𝐶can be precomputed in the offline stage. Further online effort is only required for (i) given 𝜇, solve a 𝑑× 𝒩𝐶+ 𝑑+ 𝑑2 linear system, and, (ii) given 𝜇and 𝑥, perform linear combinations and the matrix vector product in (63) employing either precomputed quantities or coefficients from (i). 2.3 Inverse distance weighting interpolation The inverse distance weighting (IDW) method has been proposed in [Shepard(1968)] to deal with interpolation of scattered data. We follow [Witteveen and Bijl(2009), Forti and Rozza(2014), Ballarin et al.(in press, 2018)Ballarin, D’Amario, Perotto, and Rozza] for its presentation and the application of IDW to shape parametrization. As in the previous section, let {𝑥𝐶𝑘}𝒩𝑐 𝑘=1 ⊂R𝑑be a set of control points. The IDW interpolant ΠIDW(𝑓) of a scalar function 𝑓: R𝑑→R is defined as ΠIDW(𝑓)(𝑥) = 𝒩𝑐 ∑︁ 𝑘=1 𝑤𝑘(𝑥) 𝑓(𝑥𝐶𝑘) 𝑥∈Ω, (69) where the weight functions 𝑤𝑘: Ω →R, for 𝑘= 1, … , 𝒩𝑐are given by 𝑤𝑘(𝑥) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ‖𝑥−𝑥𝐶𝑘‖−𝑠 𝒩𝑐 ∑︁ 𝑗=1 ‖𝑥−𝑥𝐶𝑗‖−𝑠 if 𝑥̸= 𝑥𝐶𝑘, 1 if 𝑥= 𝑥𝐶𝑘, 0 otherwise. (70) 𝑠is a positive integer, modelling the assumption that the influence of the 𝑘- th control point 𝑥𝐶𝑘on 𝑥diminishes with rate −𝑠as the distance between 𝑥and 𝑥𝐶𝑘increases. IDW interpolation trivially extends to vector functions 𝑓: R𝑑→R𝑑by application to each component 𝑓1, … , 𝑓𝑑, where the weight functions 𝑤𝑘: Ω →R do not depend on the specific component. In the case of IDW shape parametrization, for any given 𝜇, the deformed position of the control points {𝑥𝐶𝑘}𝒩𝑐 𝑘=1 is supposed to be known, and equal to 𝑦𝐶𝑘(𝜇) := 𝑓(𝑥𝐶𝑘) for 𝑘= 1, … , 𝒩𝑐. We remark that the analytic expression of 𝑓is not known, but only its action through {𝑥𝐶𝑘}𝒩𝑐 𝑘=1. This is indeed the minimum requirement to properly define (69). The deformation map is therefore ℳ(𝑥; 𝜇) = 𝒩𝑐 ∑︁ 𝑘=1 𝑤𝑘(𝑥) 𝑦𝐶𝑘(𝜇) ∀𝑥∈Ω,
Page 24
24 G. Rozza et al. In an offline-online separation effort, efficient deformation can be obtained by noting that the 𝜇dependent part is decoupled from the 𝑥-dependent weight function 𝑤𝑘(𝑥). Thus, for any 𝑥, weight terms can be precomputed once and for all and stored. The online cost of the evaluation of ℳ(𝑥; 𝜇) thus requires an inexpensive linear combination of 𝑥-dependent precomputed quantities and 𝜇-dependent control points locations. We remark that, in contrast, the RBF approach (even though still based on interpolation) required a further solution of linear system of size 𝑑× 𝒩𝐶+ 𝑑+ 𝑑2. Application in the context of fluid-structure interaction (FSI) problems between a wing (structure) and surrounding air (fluid) is shown in Figure 6. The IDW deformation of the fluid mesh resulting from a vertical displacement of the tip of the wing is depicted; the structural mesh is omitted from the picture. We refer to [Ballarin et al.(in press, 2018)Ballarin, D’Amario, Perotto, and Rozza] for more details. Fig. 6: Deformation of the fluid mesh of a fluid-structure interaction problem by IDW. 3 Beyond affinity assumptions: parametric interpolation We describe here several options to deal with cases when an exact affine decom- position of the discretized differential operators, right hand sides or outputs of interest is not existing. The section begins with a brief overview concerning the description of general non-affine problems 3.1 and later we describe the so-called empirical interpolation method (EIM) family of algorithms. This methodology be- comes particularly useful to obtain an efficient offline-online splitting also in cases
Page 25
Basic Ideas and Tools for MOR of Parametric PDEs 25 with nonlinearities and non-affine parametrization. We provide a full description of the different alternatives, starting from its standard continuous version (EIM), and presenting also its discrete (DEIM) and matrix (M-DEIM) variants. The methodologies are tested for both non-affine and non-linear problems. In section 3.2 we explain in detail the basics of the empirical interpolation method. In section 3.3 we introduce the discrete variant of the empirical interpolation method at both matrix and vector level and we mention further options to obtain an approximate affine expansion. In section 3.5 we present two examples using the EIM (section 3.5.1) and the M-DEIM algorithm to deal with both non-affinity and non-linearity (section 3.5.2). 3.1 Non-affine problems As already discussed in 1.3.4, the existence of an affine decomposition of the linear and bilinear forms of the considered problem is crucial in order to obtain a computationally efficient framework, see (26)-(28). This assumption fails to be true in several situations. Such situations occur for example in case of problems with non-affine parametric-dependency, in cases with non-linear differential operators and in cases dealing with the non-affine geometrical parametrizations introduced in section 2. In fact, in these situations, the differential operators or the right-hand sides or the outputs of interest cannot be directly written using an exact affine decom- position and we have therefore to rely on an approximate affine decomposition. The Empirical Interpolation Method (EIM) is one of the key instruments to recover an approximate affine decomposition. The EIM is a general tool for the approximation of parameterized or non- linear functions by a sum of affine terms. In the expression below we report an example for a generic parameterized function 𝑓: 𝑓(𝑥; 𝜇) ≈ 𝑄 ∑︁ 𝑞=1 𝑐𝑞(𝜇)ℎ𝑞(𝑥). (71) The EIM has been firstly proposed in [Barrault et al.(2004)Barrault, Maday, Nguyen, and Pater to deal with non-affine problems in the context of RB methods and later applied to reduced order modelling in [Grepl et al.(2007)Grepl, Maday, Nguyen, and Patera]. In [Maday et al.(2008)Maday, Nguyen, Patera, and Pau] it has been extended to a general context, a slightly different variant of EIM namely Discrete Empirical In- terpolation Method has been firstly proposed in [Chaturantabut and Sorensen(2009), Chaturantabut and Sorensen(2010)]. For more details on the a posteriori error
Page 26
26 G. Rozza et al. analysis the interested reader may see [Grepl et al.(2007)Grepl, Maday, Nguyen, and Patera, Eftang et al.(2010)Eftang, Grepl, and Patera, Chen et al.(2014)Chen, Quarteroni, and Rozz while for an extension to ℎ𝑝-adaptive EIM we refer to [Eftang and Stamm(2012)]. A generalization of the EIM family of algorithms has been proposed in [Maday and Mula(2013), Chen et al.(2014)Chen, Quarteroni, and Rozza, Maday et al.(2016) while a nonintrusive EIM technique is presented in [Casenave et al.(2014)Casenave, Ern, and L and an extension with special focus on high dimensional parameter spaces is given in [Hesthaven et al.(2014)Hesthaven, Stamm, and Zhang]. 3.2 The empirical interpolation method (EIM) The EIM is a general method to approximate a parametrized function 𝑓(𝑥; 𝜇) : Ω × 𝒫EIM →R by a linear combination of 𝑄precomputed basis functions in the case where each function 𝑓𝜇:= (·; 𝜇) belongs to some Banach space 𝒳Ω. In what follows 𝜇∈𝒫EIM is the parameter vector and 𝒫EIM is the parameter space. The EIM approximation is based on an interpolation operator I𝑄that interpolates the given function 𝑓𝜇in a set of interpolation points {𝑥𝑖}𝑄 𝑖=1 ∈Ω. The interpolant function is constructed as a linear combination of hierarchically chosen basis functions {ℎ𝑖}𝑄 𝑖=1 ∈V𝐸𝐼𝑀, where V𝐸𝐼𝑀is an approximation of the function space 𝒰that contains 𝑓, i.e. V𝐸𝐼𝑀⊆𝒰. On the contrary to other interpolation methods, that usually work with generic and multi-purpose basis functions such as polynomial functions, the EIM works with problem-specific basis functions with global support and selected hierarchically. The interpolant function can be then expressed by: I𝑄[𝑓𝜇] (𝑥) = 𝑄 ∑︁ 𝑞=1 𝑐𝑞(𝜇)ℎ𝑞(𝑥), 𝑥∈Ω, 𝜇∈𝒫𝐸𝐼𝑀, (72) where 𝑐𝑞are parameter dependent coefficients. Once the basis functions ℎ𝑞(𝑥) are set, the problem of finding the coefficients 𝑐𝑞(𝜇) is solved imposing the interpolation condition, i.e.: I𝑄[𝑓𝜇] (𝑥𝑞) = 𝑄 ∑︁ 𝑞=1 𝑐𝑞(𝜇)ℎ𝑞(𝑥𝑞) = 𝑓𝜇(𝑥𝑞), 𝑞= 1, … , 𝑄. (73) The above problem can be recast in matrix form as 𝑇𝑐𝜇= 𝑓𝜇with: (𝑇)𝑖𝑗= ℎ𝑗(𝑥𝑖), (𝑐𝜇)𝑗= 𝑐𝑗(𝜇), (𝑓(𝜇))𝑗= 𝑓(𝑥𝑖; 𝜇), 𝑖, 𝑗= 1, … , 𝑄. (74) This problem can be easily solved given the fact that the basis functions ℎ𝑞(𝑥) and the interpolation points 𝑥𝑞are known and that the matrix 𝑇is invertible.
Page 27
Basic Ideas and Tools for MOR of Parametric PDEs 27 The selection of the basis functions {ℎ𝑞}𝑄 𝑞=1 and of the interpolation points {𝑥𝑞}𝑄 𝑞=1, that are defined by a linear combination of selected function realizations {𝑓𝜇𝑖}𝑄 𝑖=1, is done following a greedy approach similar to the one presented in Subsection 1.3.2 (see Algorithm 2). The procedure provides also a set of sample points {𝜇𝑞}𝑄 𝑞=1 that are required for the construction of the basis functions. Since the basis functions are defined as linear combinations of the function realizations inside the parameter space, in order to approximate the function 𝑓 with a relatively small number of basis functions ℎ𝑞, the manifold: ℳEIM = {𝑓(𝑥; 𝜇)|𝜇∈𝒫EIM}, (75) must have a small Kolmogorov N-width [Kolmogoroff(1936)]. Once a proper norm on Ω has been defined, and here we consider 𝐿𝑝(Ω)-norms for 1 ≤𝑝≤∞, the procedure starts with the selection of the first parameter sample which is computed as: 𝜇1 = arg sup 𝜇∈𝒫EIM ‖𝑓𝜇(𝑥)‖𝐿𝑝(Ω), while the first interpolation point is computed as: 𝑥1 = arg sup 𝑥∈Ω |𝑓𝜇1(𝑥)|. The first basis function and the interpolation operator at this stage are then defined as: ℎ1(𝑥) = 𝑓𝜇1(𝑥) 𝑓𝜇1(𝑥1), I1𝑓𝜇 = 𝑓(𝑥1; 𝜇)ℎ1(𝑥). At the subsequent steps, the next basis function is selected as the one that is the worse approximated by the current interpolation operator and using a similar concept the interpolation point, often referred as magic point, is the one where the interpolation error is maximized. In mathematical terms, at the step 𝑘, the sample point is selected as the one that maximizes the error between the function 𝑓and the interpolation operator computed at the previous step I𝑘−1[𝑓]: 𝜇𝑘= arg sup 𝜇∈𝒫EIM ‖𝑓𝜇(𝑥) −I𝑘−1 [𝑓𝜇] (𝑥)‖𝐿𝑝(Ω). Once the sample point has been determined, the interpolation point is selected, in a similar fashion, as the point inside the domain that maximizes the error between the function 𝑓and the interpolation operator: 𝑥𝑘= arg sup 𝑥∈Ω |𝑓𝜇𝑘(𝑥) −I𝑘−1 [𝑓𝜇𝑘] (𝑥)|.
Page 28
28 G. Rozza et al. The next basis function is defined similarly to the first one with: ℎ𝑘(𝑥) = 𝑓𝜇𝑘(𝑥) −I𝑘−1𝑓𝜇𝑘 𝑓𝜇𝑘(𝑥𝑘) −I𝑘−1𝑓𝜇𝑘 The procedure is repeated until a certain tolerance tol is reached or a maximum number of terms 𝑁𝑚𝑎𝑥are computed (see Algorithm 2). We remark that by construction the basis functions {ℎ1, … , ℎ𝑄} and the functions {𝑓𝜇1, … , 𝑓𝜇𝑄} span the same space V𝐸𝐼𝑀: V𝐸𝐼𝑀= span{ℎ1, … , ℎ𝑄} = span{𝑓𝜇1, … , 𝑓𝜇𝑄}. However, the former are preferred for the following reasons (for more details and for the mathematical proofs we refer to [Barrault et al.(2004)Barrault, Maday, Nguyen, and P – They are linearly independent, – ℎ𝑖(𝑥𝑖) = 1 for 1 ≤𝑖≤𝑄and ℎ𝑖(𝑥𝑗) = 0 for 1 ≤𝑖≤𝑗≤𝑄, – They make the interpolation matrix 𝑇of Equation 74 to be lower triangular and with diagonal elements equal to unity and therefore the matrix is invertible. The third point implies that the interpolation problem is well-posed. Algorithm 2 The EIM algorithm - Continuous Version Input: set of parameterized functions 𝑓𝜇: Ω →R, tolerance tol and maximum number of basis functions 𝑁𝑚𝑎𝑥, 𝑝order of the chosen 𝑝-norm. Output: basis functions {ℎ1, …, ℎ𝑄}, interpolation points {𝑥1, … , 𝑥𝑄}; 𝑘= 1 ; 𝜀= tol + 1; while 𝑘< 𝑁𝑚𝑎𝑥and 𝜀> tol do Pick the Sample point: 𝜇𝑘= arg sup 𝜇∈𝒫EIM ‖𝑓𝜇(𝑥) −I𝑘−1 [𝑓𝜇] (𝑥)‖𝐿𝑝(Ω); Compute the corresponding interpolation point: 𝑥𝑘= arg sup 𝑥∈Ω |𝑓𝜇𝑘(𝑥) −I𝑘−1 [𝑓𝜇𝑘] (𝑥)|; Define the next basis function: ℎ𝑘(𝑥) = 𝑓𝜇𝑘(𝑥)−I𝑘−1𝑓𝜇𝑘 𝑓𝜇𝑘(𝑥𝑘)−I𝑘−1𝑓𝜇𝑘; Compute the error level: 𝜀= ‖𝜀𝑝‖𝐿∞with 𝜀𝑝(𝜇) = ‖𝑓𝜇(𝑥) −I𝑘−1𝑓𝜇‖𝐿𝑝(Ω); 𝑘= 𝑘+ 1; end while
Page 29
Basic Ideas and Tools for MOR of Parametric PDEs 29 3.2.1 Error Analysis Dealing with interpolation procedures, the error analysis usually involves a Lebesgue constant. In particular, in case one is using the 𝐿∞(Ω)-norm the error analysis involves the computation of the Lebesgue constant Λ𝑞= sup𝑥∈Ω ∑︀𝑞 𝑖=1 |𝐿𝑖(𝑥)| being 𝐿𝑖∈V𝐸𝐼𝑀Lagrange functions that satisfies 𝐿𝑖(𝑥𝑗) = 𝛿𝑖𝑗. It can be proved that the interpolation error is bounded by the fol- lowing expression [Barrault et al.(2004)Barrault, Maday, Nguyen, and Patera]: ‖𝑓𝜇−Iq[𝑓𝜇]‖𝐿∞(Ω) ≤(1 + Λ𝑞) inf 𝑣𝑞∈V𝐸𝐼𝑀‖𝑓𝜇−𝑣𝑞‖𝐿∞(Ω). (76) An upper bound for the Lebesgue constant, that in practice has been demon- strated to be very conservative [Barrault et al.(2004)Barrault, Maday, Nguyen, and Patera] can be computed as: Λ𝑞≤2𝑞−1. For more details concerning the estimates of the interpolation error we refer to [Barrault et al.(2004)Barrault, Maday, Nguyen, and Patera, Maday et al.(2008)Maday, Nguye 3.2.2 Practical implementation of the algorithm Practically, finding the maximum of Algorithm 2 is usually not feasible and therefore the continuous version must be transformed into a computable one. This is done selecting a finite dimensional set of training points in the parameter space {𝜇𝑖}𝑁 𝑖=1 ∈𝒫train EIM ⊂𝒫EIM and in the physical domain {𝑥𝑖}𝑀 𝑖=1 ∈ Ωℎ⊂Ω. For this reason we introduce the vector 𝑓: Ωℎ× 𝒫train EIM →R𝑀which consists into a discrete representation of the function 𝑓: (𝑓𝜇)𝑖= 𝑓𝜇(𝑥𝑖), 𝑖= 1, … , 𝑀. (77) We also define the matrix 𝐻𝑄∈R𝑀×𝑄which is defined by the discrete basis func- tions 𝐻𝑄= [ℎ1, … , ℎ𝑄] and the interpolation indices vector 𝑖𝑄= (𝑖1, … , 𝑖𝑄). The discrete interpolation operator of order 𝑄for the vector function 𝑓is then defined by: I𝑄[𝑓𝜇] = 𝐻𝑄𝑎𝑓𝜇, (78) where the coefficients 𝑎𝑓𝜇are defined such that 𝑇𝑎𝑓𝜇= 𝑓𝜇, being: 𝑇𝑘𝑞= (𝐻𝑄)𝑖𝑘𝑞, 𝑘, 𝑞= 1, … , 𝑄. (79) The implementation of the algorithm is similar to the continuous version and is reported in Algorithm 3. In the algorithm we use the notation 𝐹:,𝑗to denote the
Page 30
30 G. Rozza et al. 𝑗-th column of the matrix 𝐹. Where 𝐹∈R𝑀×𝑁is a matrix containing vector representations of the function 𝑓: (𝐹)𝑖𝑗= 𝑓(𝑥𝑖; 𝜇𝑗). (80) Once the basis and the interpolation indices are defined, during the online stage it is required to make a point-wise evaluation of the 𝑓function in the points defined by the interpolation indices. Algorithm 3 The EIM algorithm - Practical Implementation Input: set of parameter samples {𝜇𝑖}𝑀 𝑖=1 ∈𝒫train EIM ⊂𝒫EIM, set of discrete points {𝑥𝑖}𝑁 𝑖=1 ∈Ωtrain, tolerance tol, maximum number of basis functions 𝑁𝑚𝑎𝑥, 𝑝 order of the chosen 𝑝-norm. Output: basis functions matrix 𝐻𝑄= {ℎ1, …, ℎ𝑄}, interpolation indices vector 𝑖𝑄= {𝑖1, … , 𝑖𝑄}; Assemble the matrix: (𝐹)𝑖𝑗= 𝑓(𝑥𝑖; 𝜇𝑗), 𝑖= 1, … , 𝑀, 𝑗= 1, … , 𝑁; 𝑘= 1, 𝜀= tol + 1; while 𝑘< 𝑁𝑚𝑎𝑥and 𝜀> tol do Pick the sample index: 𝑗𝑘= arg max 𝑗=1,…,𝑀 ‖𝐹:,𝑗−I𝑘−1[𝐹:,𝑗]‖𝐿𝑝; and compute the interpolation point index: 𝑖𝑘= arg max 𝑖=1,…,𝑁 |𝐹𝑖,𝑗𝑘−(I𝑘−1[𝐹:,𝑗𝑘])𝑖|; define the next approximation column: ℎ𝑘= 𝐹:,𝑗𝑘−I𝑘−1[𝐹:,𝑗𝑘] 𝐹𝑖𝑘,𝑗𝑘−(I𝑘−1[𝐹:,𝑗𝑘])𝑖𝑘 define the error level: 𝜀= max 𝑗=1,…,𝑀‖𝐹:,𝑗−I𝑘−1[𝐹:,𝑗]‖𝐿𝑝 𝑘= 𝑘+ 1 end while 3.3 The Discrete Empirical Interpolation Method (DEIM) The computable version of EIM is similar to the so-called Discrete Empirical In- terpolation Method (DEIM) introduced in [Chaturantabut and Sorensen(2010)]. The main difference between the EIM and the DEIM is given by the way the basis functions are computed. In the DEIM the basis functions are computed relying on a POD procedure which is performed on a set of discrete snapshots of the parametrized function {𝑓𝑖}𝑀 𝑖=1. Each snapshot 𝑓𝑖is already considered
Page 31
Basic Ideas and Tools for MOR of Parametric PDEs 31 in discrete form in a prescribed set of points {𝑥𝑖}𝑁ℎ 𝑖=1. The procedure, which is described in detail in Algorithm 4 can be summarized into the following steps: 1. Construct the DEIM basis functions using a POD procedure on a set of previously computed snapshots: 𝐻𝑀= [ℎ1, … , ℎ𝑀] = POD(𝑓(𝜇1, … , 𝜇𝑀)). (81) 2. Given a prescribed tolerance tol determine the indices 𝑖𝑄and truncate the dimension of the POD space using an iterative greedy approach (see Algorithm 4). In Algorithm 4, with the term 𝑒𝑖𝑘, we identify a vector of dimension 𝑁ℎwhere the only non-null element is equal to 1 and is located at the index 𝑖𝑘: (𝑒𝑖𝑘)𝑗= 1 for 𝑗= 𝑖𝑘, (𝑒𝑖𝑘)𝑗= 0 for 𝑗̸= 𝑖𝑘. During the online stage, when a new value of the parameter 𝜇needs to be tested, it is required to compute the function 𝑓(𝜇) only in the location identified by the indices 𝑖𝑄. Therefore, the nonlinear function needs to be evaluated only in a relatively small number of points which is usually much smaller with respect to the total number of degrees of freedom used to discretize the domain. Algorithm 4 The DEIM procedure Input: snapshots matrix 𝑆= [𝑓(𝜇1), … , 𝑓(𝜇𝑀)] , tolerance tol. Output: DEIM basis Functions 𝐻𝑄= [ℎ1, … , ℎ𝑄], Interpolation indices 𝑖𝑄= [𝑖1, … , 𝑖𝑄]. compute the DEIM modes 𝐻𝑀= [ℎ1, … , ℎ𝑀] = POD(𝑆) 𝜀= tol + 1, 𝑘= 1 𝑖1 = arg max 𝑗=1,𝑁ℎ |(ℎ1)𝑗| 𝐻𝑄= [ℎ1], 𝑖𝑄= [𝑖1], 𝑃= [𝑒𝑖1] while 𝜀> tol do 𝑘= 𝑘+ 1 Solve (𝑃𝑇𝐻𝑄)𝑐= 𝑃𝑇ℎ𝑘 𝑟= ℎ𝑘−𝐻𝑄𝑐 𝑖𝑘= arg max 𝑗=1,𝑁ℎ |(𝑟)𝑗| 𝐻𝑄= [𝐻𝑄, ℎ𝑘], 𝑃= [𝑃, 𝑒𝑖𝑘], 𝑖𝑄= [𝑖𝑄, 𝑖𝑘] end while
Page 32
32 G. Rozza et al. 3.4 Further options Apart from the EIM and the DEIM algorithm further options are avail- able. We mention here the matrix version of the DEIM algorithm (M- DEIM) [Bonomi et al.(2017)Bonomi, Manzoni, and Quarteroni] that extends the DEIM also to the case of parametrized or non-linear matrices, the gener- alized empirical interpolation method (GEIM) [Maday and Mula(2013)] and the gappy-POD [Bui-Thanh et al.(2003)Bui-Thanh, Damodaran, and Willcox, Carlberg et al.(2010)Carlberg, Bou-Mosleh, and Farhat]. The M-DEIM is used to perform model order reduction on discretized differential operators characterized by non-linearity or non-affinity with respect to the parameter vector 𝜇. The algorithm is similar to the one in Algorithm 4 with the only difference that a vectorized version of the matrices is used to describe snapshots and POD modes. In section 3.5 we will provide an example dealing with both issues. The gappy-POD generalizes the interpolation condition to the case where the number of basis functions is smaller then the number of interpolation indices, i.e. card(𝐻𝑄) < card(𝑖𝑄). In this case the interpolation condition is substituted by a least-squares regression. The GEIM replaces the EIM requirement of a point-wise interpolation condition by the statement: 𝜎𝑗(I𝑄(𝑓(𝜇))) = 𝜎𝑗(𝑓(𝜇)), 𝑗= 1, … , 𝑄, (82) where 𝜎𝑗are a set of “well-chosen” linear functionals. For more details and for con- vergence analysis of the present method we refer to [Maday et al.(2016)Maday, Mula, and Turi 3.5 Some Examples In the previous sections we have presented the empirical interpolation method family of algorithm and we have illustrated how it is possible to recover an approximate affine expansion of the discretized differential operators. In this section we show in more detail two examples on the practical application of the EIM and the M-DEIM algorithm.
Page 33
Basic Ideas and Tools for MOR of Parametric PDEs 33 3.5.1 An heat transfer problem with a parametrized non-affine dependency forcing term In this example we illustrate the application of the computable version of the EIM on a steady state heat conduction problem in a two-dimensional square domain Ω = [−1, 1]2 with a parametrized forcing term 𝑔(𝜇) and homogeneous Dirichlet boundary conditions on the boundary 𝜕Ω. The problem is described by the following equation: {︃ −𝛼𝑡Δ𝜃= 𝑔(𝜇), in Ω, 𝜃= 0, on 𝜕Ω, (83) where 𝜃is the temperature field, 𝛼𝑡is the thermal conductivity coefficient and 𝑔(𝜇) is the parametrized forcing term which is described by the following expression: 𝑔(𝑥; 𝜇) = 𝑒−2(𝑥1−𝜇1)2−2(𝑥2−𝜇2)2, (84) where 𝜇1 and 𝜇2 are the first and second components of the parameter vector and 𝑥1 and 𝑥2 are the horizontal and vertical coordinates respectively. Let 𝑉 be a Hilbert space, the weak formulation of the problem can be written as, find 𝜃∈𝑉such that: 𝑎(𝜃(𝜇), 𝑣; 𝜇) = 𝑓(𝑣; 𝜇), ∀𝑣∈𝑉, (85) where the parametrized bilinear and linear forms are expressed by: 𝑎(𝜃, 𝑣; 𝜇) = ∫︁ Ω ∇𝜃· ∇𝑣𝑑𝑥, 𝑓(𝑣; 𝜇) = ∫︁ Ω 𝑔(𝑥; 𝜇)𝑣𝑑𝑥. (86) In the above expressions, the bilinear form 𝑎(·, ·; 𝜇) : 𝑉×𝑉→R is trivially affine while for the linear form 𝑓(·; 𝜇) : 𝑉→R we have to rely on an approximate affine expansion using the empirical interpolation method. The problem is discretized using triangular linear finite elements according to the mesh reported on left side of Figure 9. In the present case it is not possible to write an exact affine decomposition of the linear form 𝑓, we rely therefore on the computable version of the empirical interpolation method of Algorithm 3 in order to recover an approximate affine expansion. The function 𝑔(𝑥; 𝜇) is parametrized with the parameter vector 𝜇= (𝜇1, 𝜇2) ∈𝒫EIM = [−1, 1]2 that describes the position of the center of the Gaussian function. The conductivity coefficient 𝛼𝑡is fixed constant and equal to
- The testing set for the implementation of the algorithm {𝜇𝑖}𝑖=𝑁𝑡𝑟𝑎𝑖𝑛∈𝒫𝑡𝑟𝑎𝑖𝑛 is defined using 𝑁𝑡𝑟𝑎𝑖𝑛= 100 and a uniform probability distribution. The set
Page 34
34 G. Rozza et al. Fig. 7: Discretized domain into which the parameterized problem is solved (left image), together with an example of the value assumed by temperature field for one particular sample point inside the parameter space (right image). of points {𝑥𝑖}𝑁ℎ 𝑖=1 ∈Ω, that is used for the idenfitication of the magic points, is chosen to be coincident with the nodes of the finite element grid reported in Figure 7. In Figure 8 we report the first four EIM basis functions for the non-linear function 𝑔and the location of the magic points identified by the EIM algorithm. In Figure 9 we report the convergence analysis of the EIM algorithm for the nonlinear function 𝑔changing the number of EIM basis functions (left plot) and the convergence analysis of the reduced order model changing the number of reduced basis functions (right plot). 0 5 10 15 20 25 30 35 N of EIM bases 10−3 10−2 10−1 100 Average of the relative error 2 5 7 10 12 15 17 20 N of basis functions 10−2 10−1 Average of the relative error Fig. 9: Convergence analsys of the numerical example. In the left plot we can see the aver- age value of the L2 relative error between the exact function 𝑔and its EIM approximation. On the right plot we report the average value of the L2 relative error between the FOM the temparature field and the ROM temperature field. The plot is for different numbers of basis functions used to approximate the temperature field and keeping constant the number of basis functions used to approximate the forcing term (𝑁= 11)
Page 35
Basic Ideas and Tools for MOR of Parametric PDEs 35 Fig. 8: Plot of the first four modes identified by the EIM algorithm (first row and left im- age in the second row) and the location of the first 35 indices 𝑖𝑄. The magic points are identified by the red elements in the right picture on the second row. 3.5.2 An example in the context of reduced order models with non-linearity and non-affine parametric dependency In this second illustrative example we show the application of the DEIM algorithm to the stationary parametrized Navier-Stokes equations. In the present case we have both non-linearity and non-affinity with respect to the input parameters. Both non-linearity and non-affinity have been tackled using the matrix version of the Discrete Empirical Interpolation Method. The computational domain is given by the unit square Ω = [0, 1]2 and the physical problem is described by the well known Navier-Stokes equations: ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ div(𝑢⊗𝑢) −div(2𝜈(𝜇)∇𝑠𝑢) = −∇𝑝, in Ω, div 𝑢= 0, in Ω, 𝑢(𝑥) = (1, 0), on Γ𝑇𝑂𝑃, 𝑢(𝑥) = 0, on Γ0. (87) The physical problem is the classical benchmark of the lid-driven cavity problem with a parametrized diffusivity constant 𝜈(𝜇). In this case the impossibility of recovering an affine decomposition of the differential operators is given by the convective term, which is by nature a nonlinear term, and by the parametrized
Page 36
36 G. Rozza et al. diffusion term. The diffusivity constant 𝜈(𝜇) has in fact been parametrized by the following non-linear function: 𝜈(𝑥; 𝜇) = 𝑒2(−2(𝑥1−𝜇1−0.5)2−2(𝑥2−𝜇2−0.5)2) 100
- 0.01, (88) the above function is a Gaussian function and the position of its center has been parametrized using the parameter vector 𝜇= (𝜇1, 𝜇2). For the particular case, the discretized algebraic version of the continuous formulation can be rewritten as: (︂ 𝐶(𝑢) + 𝐴(𝜇) 𝐵𝑇 𝐵 0 )︂(︂ 𝑢 𝑝 )︂ = (︂ 𝑓 0 )︂ . (89) The matrix 𝐴(𝜇) represents the discretized diffusion operator, the matrix 𝐶(𝑢) represents the discretized non-linear convective operator while the term 𝐵 represents the divergence operator. The term 𝐴(𝜇) is characterized by a non-affine parametric dependency while the term 𝐶(𝑢) is characterized by non-linearity with respect to the solution. The velocity and pressure fields are approximated as: 𝑢(𝜇) ≈ 𝑁𝑢 ∑︁ 𝑞=1 𝑐𝑢 𝑞(𝜇)ℎ𝑢 𝑞, 𝑝(𝜇) ≈ 𝑁𝑝 ∑︁ 𝑞=1 𝑐𝑝 𝑞(𝜇)ℎ𝑝 𝑞, (90) and, in order to achieve an efficient offline-online splitting, the discretized op- erators are approximated by the matrix version of the DEIM algorithm and expressed as: 𝐴(𝜇) ≈ 𝑁𝐴 ∑︁ 𝑞=1 𝑐𝐴 𝑞(𝜇)ℎ𝐴 𝑞, 𝐶(𝑢) ≈ 𝑁𝐶 ∑︁ 𝑞=1 𝑐𝐶 𝑞(𝑐𝑢)ℎ𝐶 𝑞. (91) The problem is discretized using the finite volume method and a staggered cartesian grid made of 20 × 20 cell-centered finite volume elements. The DEIM algorithm has been implemented using 100 samples chosen randomly inside the training space 𝒫train EIM ∈[−0.5, 0.5]2. The magic points necessary for the implementation of the DEIM algorithm are chosen to be coincident to the cell centers of the discretized problem. The basis functions ℎ𝐴 𝑞and ℎ𝐴 𝐶are obtained using the DEIM algorithm applied on the vectorized version of the discretized differential operator snapshots computed during the training stage 𝑆𝐴= [vec(𝐴1), … , vec(𝐴𝑀)] and 𝑆𝐶= [vec(𝐶1), … , vec(𝐶𝑀)]. The snapshot matrices 𝑆𝐴and 𝑆𝐶contain in fact the discretized differential operators in vector form obtained for the different samples of the training set. In Figure 10 we report the comparison of the Full Order Model fields and the reduced order model ones; the comparison is depicted for a parameter sample
Page 37
Basic Ideas and Tools for MOR of Parametric PDEs 37 (a) (b) (c) (d) Fig. 10: Comparison between the FOM velocity (a) and pressure (c) fields and the ROM velocity (b) and pressure (d) fields. The plots are reported for one selected sample value inside the testing set. The reduced order model solutions have been computed using 14 basis functions for the velocity space, 10 for the pressure space, 10 DEIM basis functions for the convective matrix 𝐶and 10 DEIM basis functions for the diffusion matrix 𝐵. not used to train the ROM. On the right side of Figure 11 we report the convergence analysis for the numerical example. The plots are performed testing the reduced order model on 100 additional sample values selected randomly inside the parameter space 𝒫test EIM ∈[−0.5, 0.5]2. In the plots is reported the average value over the testing space of the 𝐿2 relative error. 0 5 10 15 20 25 30 N. of modes 10−9 10−7 10−5 10−3 10−1 λi λmax Conv. Term Diff. Term Pressure Velocity 0 5 10 15 20 25 30 35 40 N. of DEIM modes 10−1 6 × 10−2 L2 Relative Error Velocity Pressure Fig. 11: Eigenvalue Decay of the POD procedure during the DEIM algorithm (left plot). The convergence analysis with respect to the number of DEIM basis function (right plot), which is computed using the average value over the testing set of the 𝐿2 relative error, has been performed keeping constant the number of basis functions used to approximate the velocity and pressure fields (𝑁𝑢= 14, 𝑁𝑝= 10) and changing the number of DEIM basis functions used to approximate the convective and diffusion terms (𝑁𝐶= 𝑁𝐴).
Page 38
38 G. Rozza et al. 4 Advanced tools: reduction in parameter spaces Often the use of the aforementioned geometrical morphing techniques in Section 2 does not tell us how many control points, i.e. geometrical parameters, are enough to conduct a proper analysis. This leads to self imposing too few parameters in order to avoid the curse of dimensionality and dealing with intractable problems. To overcome this issue there exist techniques for parameter space dimensionality reduction, both linear and nonlinear. In particular we present here the active subspaces property for linear dimensionality reduction, while in the last section we show an overview of possible nonlinear methods. These methods have to be intended as general tools, not restricted to parametrised PDEs. Moreover the nature of the parameter space can be very diverse, including both geometrical and physical parameters. They are data- driven tools working with couples of input/output data, and they can be used to enhance other model order reduction techniques. 4.1 Active subspaces property and its applications In this and the following sections we present the active subspaces (AS) prop- erty proposed by Trent Russi [Russi(2010)] and developed by Paul Constan- tine [Constantine(2015)]. In brief, active subspaces are defined as the leading eigenspaces of the second moment matrix of the function’s gradient and constitute a global sensitivity index. We present how to exploit AS to reduce the parameter space dimensionality, and use it as a powerful pre-processing tool. Moreover we show how to combine it with a model reduction methodology and present its application to a cardio- vascular problem. In particular, after identifying a lower dimensional parameter subspace, we sample it to apply further model order reduction methods. This results in improved computational efficiency. The main characteristic of AS is the fact that it uses information of both the output function of interest, and the input parameter space in order to reduce its dimensionality. The active subspaces have been successfully employed in many engineering fields. We cite, among others, applications in magnetohydrodynamics power generation model in [Glaws et al.(2017)Glaws, Constantine, Shadid, and Wildey], in naval engineering for the computation of the total drag resistance with both geo- metrical and physical parameters in [Tezzele et al.(2018c)Tezzele, Salmoiraghi, Mola, and Roz Demo et al.(2018b)Demo, Tezzele, Mola, and Rozza], and for constrained shape
Page 39
Basic Ideas and Tools for MOR of Parametric PDEs 39 optimization [Lukaczyk et al.(2014)Lukaczyk, Constantine, Palacios, and Alonso] using the concept of shared active subspaces in [Tezzele et al.(2018b)Tezzele, Demo, Gadalla, Mo There are also applications to turbomachinery in [Bahamonde et al.(2017)Bahamonde, Pini, De to uncertainty quantification in the numerical simulation of a scramjet in [Constantine et al.(2015)Constantine, Emory, Larsson, and Iaccarino], and to accelerate Markov chain Monte Carlo in [Constantine et al.(2016)Constantine, Kent, and Bui-T Extension of active subspace discovery for time-dependent processes and applica- tion to a lithium ion battery model can be found in [Constantine and Doostan(2017)]. A multifidelity approach to reduce the cost of performing dimension reduc- tion through the computation of the active subspace matrix is presented in [Lam et al.(2018)Lam, Zahm, Marzouk, and Willcox]. In [Eriksson et al.(2018)Eriksson, Don they exploit AS for Bayesian optimization, while the coupling with reduced order methods can be found in [Demo et al.(2019)Demo, Tezzele, and Rozza], for a non-intrusive data-driven approach, and the coupling with POD-Galerkin methods for biomedical engineering will be presented in Section 4.4 follow- ing [Tezzele et al.(2018a)Tezzele, Ballarin, and Rozza]. 4.2 Active subspaces definition Given a parametric scalar function 𝑓(𝜇) : R𝑝→R, where 𝑝is the number of parameters, representing the output of interest, and given a probability density function 𝜌: R𝑝→R+ that represents uncertainty in the model inputs, active subspaces are low dimensional subspaces of the input space where 𝑓varies the most on average. It is a property of the of the pair (𝑓, 𝜌) (see [Constantine(2015)]). In order to uncover AS we exploit the gradients of the function with respect to the input parameters, so it can be viewed as a derivative-based sensitivity analysis that unveils low dimensional parametrization of 𝑓using some linear combinations of the original parameters. Roughly speaking, after a rescaling of the input parameter space to the hypercube [−1, 1]𝑝, we rotate it until the lower rank approximation of the output of interest is discovered, that means a preferred direction in the input space is identified. Then we can project all the data onto the orthogonal space of this preferred direction and we can construct a surrogate model on this low dimensional space. Let us add some hypotheses to 𝑓in order to proper construct the matrix we will use to find the active subspaces: let 𝑓be continuous and differentiable with square-integrable partial derivatives in the support of 𝜌. We define the so-called uncentered covariance matrix C of the gradients of 𝑓, as the matrix whose elements are the average products of partial derivatives of the map 𝑓,
Page 40
40 G. Rozza et al. that is C = E [∇𝜇𝑓∇𝜇𝑓𝑇] = ∫︁ (∇𝜇𝑓)(∇𝜇𝑓)𝑇𝜌𝑑𝜇, (92) where E is the expected value, and ∇𝜇𝑓= ∇𝑓(𝜇) = [︁ 𝜕𝑓 𝜕𝜇1 , … , 𝜕𝑓 𝜕𝜇𝑝 ]︁𝑇 is the column vector of partial derivatives of 𝑓. This matrix is symmetric so it has a real eigenvalue decomposition: C = WΛW𝑇, (93) where W ∈R𝑝×𝑝is the orthogonal matrix of eigenvectors, and Λ is the diagonal matrix of non-negative eigenvalues arranged in descending order. The eigenpairs of the uncentered covariance matrix define the active subspaces of the pair (𝑓, 𝜌). Moreover Lemma 2.1 in [Constantine et al.(2014)Constantine, Dow, and Wang] states that the eigenpairs are functionals of 𝑓(𝜇) and it holds 𝜆𝑖= w𝑇 𝑖Cw𝑖= ∫︁ (∇𝜇𝑓𝑇w𝑖)2𝜌𝑑𝜇, (94) that means that the 𝑖-th eigenvalue is the average squared directional derivative of 𝑓along the eigenvector w𝑖. Alternatively we can say that the eigenvalues represent the magnitude of the variance of ∇𝜇𝑓along their eigenvectors orientation. So small values of the eigenvalues correspond to small perturbation of 𝑓along the corresponding eigenvectors. It also follows that large gaps between eigenvalues indicate directions where 𝑓changes the most on average. Since we seek for a lower dimensional space of dimension 𝑀< 𝑝where the target function has exactly this property, we define the active subspace of dimension 𝑀as the span of the first 𝑀eigenvectors (they correspond to the most energetic eigenvalues before a gap). Let us partition Λ and W as Λ = [︂ Λ1 Λ2 ]︂ , W = [W1 W2] , (95) where Λ1 = diag(𝜆1, … , 𝜆𝑀), and W1 contains the first 𝑀eigenvectors. We can use W1 to project the original parameters to the active subspace obtaining the reduced parameters, that is the input space is geometrically transformed and aligned with W1, in order to retain only the directions where the function variability is high. We call active variable 𝜇𝑀the range of W𝑇 1 , and inactive variable 𝜂the range of W𝑇 2 : 𝜇𝑀= W𝑇 1 𝜇∈R𝑀, 𝜂= W𝑇 2 𝜇∈R𝑝−𝑀. (96) We can thus express any point in the parameter space 𝜇∈R𝑝in terms of 𝜇𝑀 and 𝜂as 𝜇= WW𝑇𝜇= W1W𝑇 1 𝜇+ W2W𝑇 2 𝜇= W1𝜇𝑀+ W2𝜂. (97)
Page 41
Basic Ideas and Tools for MOR of Parametric PDEs 41 The lower-dimension approximation, or surrogate quantity of interest, 𝑔: R𝑀→R of the target function 𝑓, is a function of only the active variable 𝜇𝑀as 𝑓(𝜇) ≈𝑔(W𝑇 1 𝜇) = 𝑔(𝜇𝑀). (98) Such 𝑔is called ridge function (see [Pinkus(2015)]) and, as we can infer from this section, it is constant along the span of W2. From a practical point of view Eq. (92) is estimated through Monte Carlo method. We draw 𝑁train independent samples 𝜇(𝑖) according to the measure 𝜌 and we approximate C ≈^C = 1 𝑁train 𝑁train ∑︁ 𝑖=1 ∇𝜇𝑓𝑖∇𝜇𝑓𝑇 𝑖= ^ W ^Λ ^ W𝑇, (99) where ∇𝜇𝑓𝑖= ∇𝜇𝑓(𝜇(𝑖)). In [Constantine(2015)] they provide an heuristic formula for the number of samples 𝑁train needed to properly estimate the first 𝑘 eigenvalues, that is 𝑁train = 𝛼𝑘ln(𝑝), (100) where 𝛼usually is between 2 and 10. Moreover they prove that for sufficiently large 𝑁train the error 𝜀committed in the approximation of the active subspace of dimension 𝑛is bounded from above by 𝜀= dist(rank(W1), rank( ^ W1)) ≤ 4𝜆1𝛿 𝜆𝑛−𝜆𝑛+1 , (101) where 𝛿is a positive scalar bounded from above by 𝜆𝑛−𝜆𝑛+1 5𝜆1 . Here we can clearly see how the gap between two eigenvalues is important in order to properly approximate 𝑓exploiting AS. 4.3 Some examples In this section we present two simple examples with the computation of the active subspaces using analytical gradients. To highlight the possibility that the presence of an active subspace is not alway guaranteed we also show an example in this direction. We choose for both the cases a tridimensional input parameter space without loss of generality. In order to identify the low-dimensional structure of the function of interest we use the sufficient summary plots, developed in [Cook(2009)]. In our cases, they are scatter plots of 𝑓(𝜇) against the active variable 𝜇𝑀. The presence of an active subspace is not always guaranteed. For example every target function that has a radial symmetry has not a lower dimensional
Page 42
42 G. Rozza et al. representation in terms of active variable. This is due to the fact that there is no rotation of the input parameter space that aligns it along a preferred direction since all of them are equally important. Let us consider for example the function 𝑓(𝜇) = 1 2𝜇𝑇𝜇representing an 𝑛-dimensional elliptic paraboloid, where the parameter 𝜇is a column vector in [−1, 1]3. In this case we have the exact derivatives, in fact ∇𝜇𝑓= 𝜇and we do not have to approximate them. If we draw 1000 samples and we apply the procedure to find an active subspace and we plot the sufficient summary plot in one dimension, as in Figure 12, we clearly see how it is unable to find the active variable along which 𝑓varies the most on average. In fact there is not a significant gap between the eigenvalues, since we have that C = 1 3Id. Moreover the projection of the data onto the inactive subspace suggest us the presence of an 𝑛-dimensional elliptic paraboloid. Fig. 12: Example of an output function with a radial symmetry. On the left the exact eigen- values of the uncentered covariance matrix. On the right the sufficient summary plot in one dimension (𝑓(𝜇) against 𝜇𝑀= W𝑇 1 𝜇) shows how the projection of the data along the inactive directions does not unveil a lower dimensional structure for 𝑓. Let us consider now another quadratic function in 3 variables. We define the output of interest 𝑓as 𝑓(𝜇) = 1 2𝜇𝑇A𝜇, (102) where 𝜇∈[−1, 1]3, and A is symmetric positive definite with a major gap between the first and the second eigenvalue. With this form we can compute the exact gradients as ∇𝜇𝑓(𝜇) = A𝜇and, taking 𝜌as a uniform density function,
Page 43
Basic Ideas and Tools for MOR of Parametric PDEs 43 compute C as C = A (︂∫︁ 𝜇𝜇𝑇𝜌𝑑𝜇 )︂ A𝑇= 1 3A2. (103) So the squared eigenvalues of A are the eigenvalues of C. Since, by definition, A has a significant gap between the first and the second eigenvalues we can easily find an active subspace of dimension one. In Figure 13 we show the sufficient summary plot of 𝑓with respect to its active variable. A clear univariate behavior is present, as expected, so we can easily construct 𝑔, for instance taking a quadratic unidimensional function. We can also see the associated eigenvalues of the uncentered covariance matrix. Fig. 13: Example of a quadratic function with an active subspace of dimension one. On the left the exact eigenvalues of the uncentered covariance matrix. On the right the sufficient summary plot in one dimension (𝑓(𝜇) against 𝜇𝑀= W𝑇 1 𝜇) shows how the projection of the data along the inactive directions unveils a univariate structure for 𝑓. 4.4 Active subspaces as pre-processing tool to enhance model reduction The presence of an active subspace for an output of interest, derived from the solution of a parametric PDE, can be exploited for further model order reduction. Thus, in this context, AS can be seen as a powerful pre-processing technique to both reduce the parameter space dimensionality and boost the performance of other model order reduction methods. In [Tezzele et al.(2018a)Tezzele, Ballarin, and Rozza] the active subspace for a relative pressure drop in a stenosed carotid artery is used as a reduced
Page 44
44 G. Rozza et al. sampling space to improve the reconstruction of the output manifold. We used as parameters the displacement of a selection radial basis functions (RBF) control points to simulate the occlusion of the carotid artery after the bifurcation. For a review of RBF interpolation technique see Section 2.2. In Figure 5 two different views of the same carotid and the control points highlighted with green dots. The target function was a relative pressure drop between the two branches computed solving a stationary Navier-Stokes problem. After the identification of the active subspace we exploit it by sampling the original full parameter space along the active subspace. These sampled parameters were used, in the offline phase, to construct the snapshots matrix for the training of a ROM. This leads to better approximation properties for a given number of snapshots with respect to usual sampling techniques. The natural construction of the uncentered covariance matrix, which uses information from both the inputs and the outputs is the reason of such improvements. The same idea has been coupled also with non-intrusive model order reduction techniques, such as proper orthogonal decomposition with interpolation (PODI), in [Tezzele et al.(2019)Tezzele, Demo, and Rozza], while for the reconstruction of modal coefficient using PODI with AS for low computational budget we suggest [Demo et al.(2019)Demo, Tezzele, and Rozza]. 4.5 About nonlinear dimensionality reduction There are plenty of other techniques that reduce the dimensionality of a given dataset. They do not exploit simultaneously the structure of the output function and the input parameter space like AS, they just express the datavectors we want to reduce in a reduced space embedded in the original one. For a comprehensive overview see [Lee and Verleysen(2007)] and [Van Der Maaten et al.(2009)Van Der Maaten, Postma, and Van den Herik]. The main assumption is that the dataset at hand has an intrinsic dimension- ality, which is lower than the full space where they belong. This means that the data are lying on or near a manifold with dimensionality 𝑑embedded in a greater space of dimension 𝐷. If we approximate this manifold with a linear subspace we use a linear dimensionality reduction technique, oth- erwise assuming the data lie on a curved manifold we can achieve better results using a nonlinear method. Unfortunately in general neither the char- acteristics of the manifold, nor the intrinsic dimensionality are known, so the dimensionality reduction problem is ill-posed. There are several algo- rithms to detect the intrinsic dimensionality of a dataset, we suggest the review in [Camastra(2003)]. Among all we cite two of the most popular techniques, which
Page 45
Basic Ideas and Tools for MOR of Parametric PDEs 45 are locally linear embedding (LLE) presented in [Roweis and Saul(2000)], and Isomap [Tenenbaum et al.(2000)Tenenbaum, De Silva, and Langford]. Exten- sions for the two methods can be found in [Bengio et al.(2004)Bengio, Paiement, Vincent, Delall LLE seeks to preserve local properties of the high dimensional data in the embedded space, and it is able to detect non-convex manifolds. In particular it preserves local reconstruction weights of the neighborhood graph, that is LLE fits a hyperplane through each data point and its nearest neighbors. Some applica- tions can be found in [González et al.(2016)González, Cueto, and Chinesta] for biomedical engineering, or in [Ibanez et al.(2018)Ibanez, Abisset-Chavanne, Aguado, Gonzalez, for computational mechanics. Isomap instead seeks to preserve geodesic (or curvilinear) distances between the high dimensional data points and the lower dimensional embedded ones. Its topological stability has been investigated in [Balasubramanian and Schwartz(2002)], while it has been used in for micromotility reconstruction in [Arroyo et al.(2012)Arroyo, Heltai, M Other approaches include for example a manifold walking algorithm that has been proposed in [Meng et al.(2015)Meng, Breitkopf, Raghavan, Mauvoisin, Bartier, and Herno and in [Meng et al.(2018)Meng, Breitkopf, Le Quilliec, Raghavan, and Villon]. 5 Conclusion and Outlook This introductory chapter provided the means to understand projection based MOR methods in section 1. Various techniques allowing the parametrization of complicated geometries are provided in section 2. Since many geometries of interest introduce non-linearities or non-affine parameter dependency an intermediate step such as the empirical interpolation method is often applied. The basics were presented in section 3 and will be used further in the following chapters. The reduction in parameter space becomes necessary if high-dimensional parameter space are considered. Active subspaces (see section 4) provide a mean to tackle the curse of dimensionality. Each chapter of the handbook gives an in-depth technical details upon a particular topic of interest. This includes common MOR methods, several application areas of interest and a survey of current software frameworks for model reduction. Whenever a method does not rely only on the PDE-based functional analysis setting introduced in this chapter, corresponding requirements will be mentoined within each technical chapter. Acknowledgment: We are grateful to EU-COST European Union Cooperation in Science and Technology, section EU-MORNET Model Reduction Network,
Page 46
46 REFERENCES TD 1307 for pushing us into this initiative. This work is supported by European Union Funding for Research and Innovation — Horizon 2020 Program — in the framework of European Research Council Executive Agency: H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” P.I. Gianluigi Rozza. References [Adams(1975)] Adams, R. A. (1975): Sobolev Spaces, Academic Press, New York. [Ainsworth and Oden(1997)] Ainsworth, M. and J. Oden (1997): “A posteriori error estimation in finite element analysis,” Computer Methods in Applied Mechanics and Engineering, 142, 1 – 88. [Almroth et al.(1978)Almroth, Stern, and Brogan] Almroth, B., P. Stern, and F. Brogan (1978): “Automatic choice of global shape functions in structural analysis,” AIAA J., 16, 525–528. [Arnold et al.(2002)Arnold, Brezzi, Cockburn, and Marini] Arnold, D., F. Brezzi, B. Cockburn, and L. Marini (2002): “Unified analysis of Discontinuous Galerkin methods for elliptic problems,” SIAM Journal on Numerical Analysis, 39, 1749–1779. [Arroyo et al.(2012)Arroyo, Heltai, Millán, and DeSimone] Arroyo, M., L. Heltai, D. Millán, and A. DeSimone (2012): “Reverse engineer- ing the euglenoid movement,” Proceedings of the National Academy of Sciences, 109, 17874–17879. [Babuska(1970/71)] Babuska, I. (1970/71): “Error-bounds for finite element method,” Numerische Mathematik, 16, 322–333. [Bahamonde et al.(2017)Bahamonde, Pini, De Servi, and Colonna] Bahamonde, S., M. Pini, C. De Servi, and P. Colonna (2017): “Active sub- spaces for the optimal meanline design of unconventional turbomachinery,” Applied Thermal Engineering, 127, 1108–1118. [Balasubramanian and Schwartz(2002)] Balasubramanian, M. and E. L. Schwartz (2002): “The isomap algorithm and topological stability,” Science, 295, 7–7. [Ballarin et al.(in press, 2018)Ballarin, D’Amario, Perotto, and Rozza] Ballarin, F., A. D’Amario, S. Perotto, and G. Rozza (in press, 2018): “A POD-Selective Inverse Distance Weighting method for fast parametrized shape morphing,” Int. J. Num. Meth. Eng. [Barrault et al.(2004)Barrault, Maday, Nguyen, and Patera] Barrault, M., Y. Maday, N. C. Nguyen, and A. T. Patera (2004): “An ‘empirical
Page 47
REFERENCES 47 interpolation’ method: application to efficient reduced-basis discretization of partial differential equations,” Comptes Rendus Mathematique, 339, 667–672. [Beckert and Wendland(2001)] Beckert, A. and H. Wendland (2001): “Multivari- ate interpolation for fluid-structure-interaction problems using radial basis functions,” Aerospace Science and Technology, 5, 125–134. [Bengio et al.(2004)Bengio, Paiement, Vincent, Delalleau, Roux, and Ouimet] Bengio, Y., J.-f. Paiement, P. Vincent, O. Delalleau, N. L. Roux, and M. Ouimet (2004): “Out-of-sample extensions for lle, isomap, mds, eigenmaps, and spectral clustering,” in Advances in neural information processing systems, 177–184. [Boffi et al.(2013)Boffi, Brezzi, and Fortin] Boffi, D., F. Brezzi, and M. Fortin (2013): Mixed Finite Element Methods and Applications, Springer-Verlag Berlin Heidelberg. [Bonomi et al.(2017)Bonomi, Manzoni, and Quarteroni] Bonomi, D., A. Man- zoni, and A. Quarteroni (2017): “A matrix DEIM technique for model reduction of nonlinear parametrized problems in cardiac mechanics,” Com- puter Methods in Applied Mechanics and Engineering, 324, 300–326. [Buhmann(2003)] Buhmann, M. D. (2003): Radial basis functions: theory and implementations, volume 12, Cambridge university press. [Bui-Thanh et al.(2003)Bui-Thanh, Damodaran, and Willcox] Bui-Thanh, T., M. Damodaran, and K. Willcox (2003): “Proper Orthogonal Decomposi- tion extensions for parametric applications in compressible aerodynamics,” in 21st AIAA Applied Aerodynamics Conference, American Institute of Aeronautics and Astronautics. [Camastra(2003)] Camastra, F. (2003): “Data dimensionality estimation meth- ods: a survey,” Pattern recognition, 36, 2945–2954. [Canuto et al.(2006)Canuto, Hussaini, Quarteroni, and Zhang] Canuto, C., M. Y. Hussaini, A. Quarteroni, and T. Zhang (2006): Spectral Methods: Fundamentals in Single Domains, Springer-Verlag Berlin Heidelberg. [Canuto et al.(2009)Canuto, Tonn, and Urban] Canuto, C., T. Tonn, and K. Ur- ban (2009): “A posteriori error analysis of the reduced basis method for nonaffine parametrized nonlinear pdes,” SIAM Journal on Numerical Anal- ysis, 47, 2001–2022. [Carlberg et al.(2010)Carlberg, Bou-Mosleh, and Farhat] Carlberg, K., C. Bou- Mosleh, and C. Farhat (2010): “Efficient non-linear model reduction via a least-squares Petrov-Galerkin projection and compressive tensor approx- imations,” International Journal for Numerical Methods in Engineering, 86, 155–181. [Casenave et al.(2014)Casenave, Ern, and Lelièvre] Casenave, F., A. Ern, and
Page 48
48 REFERENCES T. Lelièvre (2014): “A nonintrusive reduced basis method applied to aeroa- coustic simulations,” Advances in Computational Mathematics, 41, 961–986. [Chaturantabut and Sorensen(2009)] Chaturantabut, S. and D. C. Sorensen (2009): “Discrete empirical interpolation for nonlinear model reduction,” in Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, IEEE. [Chaturantabut and Sorensen(2010)] Chaturantabut, S. and D. C. Sorensen (2010): “Nonlinear model reduction via discrete empirical interpolation,” SIAM Journal on Scientific Computing, 32, 2737–2764. [Chen et al.(2014)Chen, Quarteroni, and Rozza] Chen, P., A. Quarteroni, and G. Rozza (2014): “A weighted empirical interpolation method: a priori convergence analysis and applications,” ESAIM: Mathematical Modelling and Numerical Analysis, 48, 943–953. [Chen(2016)] Chen, Y. (2016): “A certified natural-norm successive constraint method for parametric inf–sup lower bounds,” Applied Numerical Mathe- matics, 99, 98–108. [Chen et al.(2008)Chen, Hesthaven, Maday, and Rodríguez] Chen, Y., J. S. Hes- thaven, Y. Maday, and J. Rodríguez (2008): “A monotonic evaluation of lower bounds for inf-sup stability constants in the frame of reduced basis approximations,” Comptes Rendus Mathematique, 346, 1295–1300. [Chen et al.(2009)Chen, Hesthaven, Maday, and Rodríguez] Chen, Y., J. S. Hes- thaven, Y. Maday, and J. Rodríguez (2009): “Improved successive constraint method based a posteriori error estimate for reduced basis approximation of 2d Maxwell’s problem,” ESAIM: Mathematical Modelling and Numerical Analysis, 43, 1099–1116. [Chinesta et al.(2017)Chinesta, Huerta, Rozza, and Willcox] Chinesta, F., A. Huerta, G. Rozza, and K. Willcox (2017): “Model order reduction,” Encyclopedia of Computational Mechanics Second Edition, Elsevier, 1–36. [Ciarlet(2002)] Ciarlet, P. (2002): The Finite Element Method for Elliptic Prob- lems, volume 40, Classics in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia. [Ciarlet(2014)] Ciarlet, P. (2014): Linear and nonlinear functional analysis with applications, Society for Industrial and Applied Mathematics, Philadelphia. [Constantine(2015)] Constantine, P. G. (2015): Active subspaces: Emerging ideas for dimension reduction in parameter studies, volume 2, SIAM. [Constantine and Doostan(2017)] Constantine, P. G. and A. Doostan (2017): “Time-dependent global sensitivity analysis with active subspaces for a lithium ion battery model,” Statistical Analysis and Data Mining: The ASA Data Science Journal, 10, 243–262. [Constantine et al.(2014)Constantine, Dow, and Wang] Constantine, P. G.,
Page 49
REFERENCES 49 E. Dow, and Q. Wang (2014): “Active subspace methods in theory and practice: applications to kriging surfaces,” SIAM Journal on Scientific Computing, 36, A1500–A1524. [Constantine et al.(2015)Constantine, Emory, Larsson, and Iaccarino] Constantine, P. G., M. Emory, J. Larsson, and G. Iaccarino (2015): “Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet,” Journal of Computational Physics, 302, 1–20. [Constantine et al.(2016)Constantine, Kent, and Bui-Thanh] Constantine, P. G., C. Kent, and T. Bui-Thanh (2016): “Accelerating Markov chain Monte Carlo with active subspaces,” SIAM Journal on Scientific Computing, 38, A2779–A2805. [Cook(2009)] Cook, R. D. (2009): Regression graphics: ideas for studying regres- sions through graphics, volume 482, John Wiley & Sons. [Demo et al.(2018a)Demo, Tezzele, Gustin, Lavini, and Rozza] Demo, N., M. Tezzele, G. Gustin, G. Lavini, and G. Rozza (2018a): “Shape optimization by means of proper orthogonal decomposition and dynamic mode decomposition,” in Technology and Science for the Ships of the Future: Proceedings of NAV 2018: 19th International Conference on Ship & Maritime Research, IOS Press, 212–219. [Demo et al.(2018b)Demo, Tezzele, Mola, and Rozza] Demo, N., M. Tezzele, A. Mola, and G. Rozza (2018b): “An efficient shape parametrisation by free- form deformation enhanced by active subspace for hull hydrodynamic ship design problems in open source environment,” in The 28th International Ocean and Polar Engineering Conference. [Demo et al.(2019)Demo, Tezzele, and Rozza] Demo, N., M. Tezzele, and G. Rozza (2019): “A non-intrusive approach for proper orthogonal de- composition modal coefficients reconstruction through active subspaces,” Comptes Rendus de l’Académie des Sciences DataBEST 2019 Special Issue. [Deparis et al.(2014)Deparis, Forti, and Quarteroni] Deparis, S., D. Forti, and A. Quarteroni (2014): “A rescaled localized radial basis function inter- polation on non-Cartesian and nonconforming grids,” SIAM Journal on Scientific Computing, 36, A2745–A2762. [Deparis and Rozza(2009)] Deparis, S. and G. Rozza (2009): “Reduced basis method for multi-parameter-dependent steady navier–stokes equations: Applications to natural convection in a cavity,” Journal of Computational Physics, 228, 4359–4378. [Duchon(1977)] Duchon, J. (1977): “Splines minimizing rotation-invariant semi- norms in Sobolev spaces,” Constructive theory of functions of several variables, 85–100.
Page 50
50 REFERENCES [Eckart and Young(1936)] Eckart, C. and G. Young (1936): “The approximation of one matrix by another of lower rank. psychometrika,” Psychometrika, 1, 211–218. [Eftang et al.(2010)Eftang, Grepl, and Patera] Eftang, J. L., M. A. Grepl, and A. T. Patera (2010): “A posteriori error bounds for the empirical interpo- lation method,” Comptes Rendus Mathematique, 348, 575–579. [Eftang and Stamm(2012)] Eftang, J. L. and B. Stamm (2012): “Parameter multi-domain ‘hp’ empirical interpolation,” International Journal for Nu- merical Methods in Engineering, 90, 412–428. [Eriksson et al.(2018)Eriksson, Dong, Lee, Bindel, and Wilson] Eriksson, D., K. Dong, E. H. Lee, D. Bindel, and A. G. Wilson (2018): “Scaling Gaussian process regression with derivatives,” arXiv preprint arXiv:1810.12283. [Evans(1998)] Evans, L. (1998): Partial Differential Equations, American Math- ematical Society. [Eymard et al.(2000)Eymard, Gallouët, and Herbin] Eymard, R., T. R. Gal- louët, and R. Herbin (2000): “The finite volume method,” in Handbook of Numerical Analysis, volume 7, Editors: P.G. Ciarlet and J.L. Lions, 713–1020. [Forti and Rozza(2014)] Forti, D. and G. Rozza (2014): “Efficient geometrical parametrisation techniques of interfaces for reduced-order modelling: ap- plication to fluid–structure interaction coupling problems,” International Journal of Computational Fluid Dynamics, 28, 158–169. [Glaws et al.(2017)Glaws, Constantine, Shadid, and Wildey] Glaws, A., P. G. Constantine, J. N. Shadid, and T. M. Wildey (2017): “Dimension re- duction in magnetohydrodynamics power generation models: Dimensional analysis and active subspaces,” Statistical Analysis and Data Mining: The ASA Data Science Journal, 10, 312–325. [González et al.(2016)González, Cueto, and Chinesta] González, D., E. Cueto, and F. Chinesta (2016): “Computational patient avatars for surgery plan- ning,” Annals of biomedical engineering, 44, 35–45. [Graetsch and Bathe(2005)] Graetsch, T. and K.-J. Bathe (2005): “A posteriori error estimation techniques in practical finite element analysis,” Computers & Structures, 83, 235 – 265. [Grepl et al.(2007)Grepl, Maday, Nguyen, and Patera] Grepl, M. A., Y. Maday, N. C. Nguyen, and A. T. Patera (2007): “Efficient reduced-basis treat- ment of nonaffine and nonlinear partial differential equations,” ESAIM: Mathematical Modelling and Numerical Analysis, 41, 575–605. [Hess et al.(2015)Hess, Grundel, and Benner] Hess, M., S. Grundel, and P. Ben- ner (2015): “Estimating the inf-sup constant in reduced basis methods for
Page 51
REFERENCES 51 time-harmonic Maxwell’s equations,” IEEE Transactions on Microwave Theory and Techniques, 63, 3549–3557. [Hesthaven et al.(2012)Hesthaven, Stamm, and Zhang] Hesthaven, J., B. Stamm, and S. Zhang (2012): “Certified reduced basis method for the electric field integral equation,” SIAM Journal on Scientific Computing, 34, A1777–A1799. [Hesthaven et al.(2014)Hesthaven, Stamm, and Zhang] Hesthaven, J. S., B. Stamm, and S. Zhang (2014): “Efficient greedy algorithms for high- dimensional parameter spaces with applications to empirical interpolation and reduced basis methods,” ESAIM: Mathematical Modelling and Numerical Analysis, 48, 259–283. [Huynh et al.(2007)Huynh, Rozza, Sen, and Patera] Huynh, D., G. Rozza, S. Sen, and A. Patera (2007): “A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants,” Comptes Rendus Mathematique, 345, 473–478. [Huynh et al.(2010)Huynh, Knezevic, Chen, Hesthaven, and Patera] Huynh, D. B. P., D. J. Knezevic, Y. Chen, J. S. Hesthaven, and A. T. Patera (2010): “A natural-norm successive constraint method for inf-sup lower bounds,” Computer Methods in Applied Mechanics and Engineering, 199, 1963–1975. [Iapichino et al.(2017)Iapichino, Ulbrich, and Volkwein] Iapichino, L., S. Ul- brich, and S. Volkwein (2017): “Multiobjective PDE-constrained opti- mization using the reduced-basis method,” Advances in Computational Mathematics, 43, 945–972. [Ibanez et al.(2018)Ibanez, Abisset-Chavanne, Aguado, Gonzalez, Cueto, and Chinesta] Ibanez, R., E. Abisset-Chavanne, J. V. Aguado, D. Gonzalez, E. Cueto, and F. Chinesta (2018): “A manifold learning approach to data-driven computational elasticity and inelasticity,” Archives of Computational Methods in Engineering, 25, 47–57. [Kolmogoroff(1936)] Kolmogoroff, A. (1936): “Uber die beste annaherung von funktionen einer gegebenen funktionenklasse,” The Annals of Mathematics, 37, 107, URL https://doi.org/10.2307/1968691. [Lam et al.(2018)Lam, Zahm, Marzouk, and Willcox] Lam, R., O. Zahm, Y. Marzouk, and K. Willcox (2018): “Multifidelity dimension reduction via active subspaces,” arXiv preprint arXiv:1809.05567. [Lassila and Rozza(2010)] Lassila, T. and G. Rozza (2010): “Parametric free- form shape design with PDE models and reduced basis method,” Computer Methods in Applied Mechanics and Engineering, 199, 1583–1592. [Lee and Verleysen(2007)] Lee, J. A. and M. Verleysen (2007): Nonlinear dimen- sionality reduction, Springer Science & Business Media. [Lombardi et al.(2012)Lombardi, Parolini, Quarteroni, and Rozza] Lombardi,
Page 52
52 REFERENCES M., N. Parolini, A. Quarteroni, and G. Rozza (2012): “Numerical simulation of sailing boats: Dynamics, FSI, and shape optimization,” in Variational Analysis and Aerospace Engineering: Mathematical Challenges for Aerospace Design, Springer, 339. [Lukaczyk et al.(2014)Lukaczyk, Constantine, Palacios, and Alonso] Lukaczyk, T. W., P. Constantine, F. Palacios, and J. J. Alonso (2014): “Active subspaces for shape optimization,” in 10th AIAA multidisciplinary design optimization conference, 1171. [Maday and Mula(2013)] Maday, Y. and O. Mula (2013): “A generalized empir- ical interpolation method: Application of reduced basis techniques to data assimilation,” in Analysis and Numerics of Partial Differential Equations, Springer Milan, 221–235. [Maday et al.(2016)Maday, Mula, and Turinici] Maday, Y., O. Mula, and G. Turinici (2016): “Convergence analysis of the generalized empirical interpolation method,” SIAM Journal on Numerical Analysis, 54, 1713– 1731. [Maday et al.(2008)Maday, Nguyen, Patera, and Pau] Maday, Y., N. Nguyen, A. Patera, and S. Pau (2008): “A general multipurpose interpolation pro- cedure: the magic points,” Communications on Pure and Applied Analysis, 8, 383–404. [Manzoni(2014)] Manzoni, A. (2014): “An efficient computational framework for reduced basis approximation and a posteriori error estimation of parametrized navier–stokes flows,” ESAIM: Mathematical Modelling and Numerical Analysis, 48, 1199–1226. [Manzoni and Negri(2015)] Manzoni, A. and F. Negri (2015): “Heuristic strate- gies for the approximation of stability factors in quadratically nonlinear parametrized pdes,” Advances in Computational Mathematics, 41, 1255– 1288. [Manzoni et al.(2012)Manzoni, Quarteroni, and Rozza] Manzoni, A., A. Quar- teroni, and G. Rozza (2012): “Model reduction techniques for fast blood flow simulation in parametrized geometries,” International journal for numerical methods in biomedical engineering, 28, 604–625. [Meng et al.(2018)Meng, Breitkopf, Le Quilliec, Raghavan, and Villon] Meng, L., P. Breitkopf, G. Le Quilliec, B. Raghavan, and P. Villon (2018): “Non- linear shape-manifold learning approach: concepts, tools and applications,” Archives of Computational Methods in Engineering, 25, 1–21. [Meng et al.(2015)Meng, Breitkopf, Raghavan, Mauvoisin, Bartier, and Hernot] Meng, L., P. Breitkopf, B. Raghavan, G. Mauvoisin, O. Bartier, and X. Hernot (2015): “Identification of material properties using indentation
Page 53
REFERENCES 53 test and shape manifold learning approach,” Computer Methods in Applied Mechanics and Engineering, 297, 239–257. [Morris et al.(2008)Morris, Allen, and Rendall] Morris, A., C. Allen, and T. Ren- dall (2008): “CFD-based optimization of aerofoils using radial basis func- tions for domain element parameterization and mesh deformation,” Inter- national Journal for Numerical Methods in Fluids, 58, 827–860. [Noor(1982)] Noor, A. (1982): “On making large nonlinear problems small,” Comput. Meth. Appl. Mech. Engrg., 34, 955–985. [Pinkus(2015)] Pinkus, A. (2015): Ridge functions, volume 205, Cambridge Uni- versity Press. [Prud’Homme et al.(2002)Prud’Homme, Rovas, Veroy, Machiels, Maday, Patera, and Turinici] Prud’Homme, C., D. V. Rovas, K. Veroy, L. Machiels, Y. Maday, A. T. Patera, and G. Turinici (2002): “Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods,” Journal of Fluids Engineering, 124, 70–80. [PyGeM(2017)] PyGeM (2017): “Python Geometrical Morphing,” URL https: //github.com/mathLab/PyGeM. [Quarteroni(2017)] Quarteroni, A. (2017): Numerical Models for Differential Problems, Modeling, Simulation and Applications, volume 16, Springer International Publishing. [Rebollo et al.(2017)Rebollo, Ávila, Mármol, Ballarin, and Rozza] Rebollo, T., E. Ávila, M. Mármol, F. Ballarin, and G. Rozza (2017): “On a certified Smagorinsky reduced basis turbulence model,” SIAM Journal on Numerical Analysis, 55, 3047–3067. [Renardy and Rogers(2004)] Renardy, M. and R. Rogers (2004): An Introduction to Partial Differential Equations, Springer-Verlag New York. [Roweis and Saul(2000)] Roweis, S. T. and L. K. Saul (2000): “Nonlinear dimen- sionality reduction by locally linear embedding,” science, 290, 2323–2326. [Rozza et al.(2008)Rozza, Huynh, and Patera] Rozza, G., D. Huynh, and A. T. Patera (2008): “Reduced Basis Approximation and a Posteriori Error Estimation for Affinely Parametrized Elliptic Coercive Partial Differential Equations,” Archives of Computational Methods in Engineering, 15, 229 – 275. [Rozza et al.(2009)Rozza, Huynh, Nguyen, and Patera] Rozza, G., D. B. P. Huynh, N. C. Nguyen, and A. T. Patera (2009): “Real-time reliable simula- tion of heat transfer phenomena,” in ASME -American Society of Mechan- ical Engineers - Heat Transfer Summer Conference, paper HT2009-88212, volume 3, 851–860. [Rozza et al.(2013)Rozza, Koshakji, and Quarteroni] Rozza, G., A. Koshakji, and A. Quarteroni (2013): “Free Form Deformation techniques applied
Page 54
54 REFERENCES to 3D shape optimization problems,” Communications in Applied and Industrial Mathematics, 4, 1–26. [Rozza et al.(2018)Rozza, Malik, Demo, Tezzele, Girfoglio, Stabile, and Mola] Rozza, G., M. H. Malik, N. Demo, M. Tezzele, M. Girfoglio, G. Stabile, and A. Mola (2018): “Advances in Reduced Order Methods for Parametric Industrial Problems in Computational Fluid Dynamics,” Glasgow, UK: ECCOMAS ECCM - ECFD Conference Proceedings. [Rudin(1976)] Rudin, W. (1976): Principles of Mathematical Analysis, Interna- tional Series in Pure and Applied Mathematics. McGraw-Hill. [Russi(2010)] Russi, T. M. (2010): Uncertainty quantification with experimental data and complex system models, Ph.D. thesis, UC Berkeley. [Salmoiraghi et al.(2016)Salmoiraghi, Ballarin, Corsi, Mola, Tezzele, and Rozza] Salmoiraghi, F., F. Ballarin, G. Corsi, A. Mola, M. Tezzele, and G. Rozza (2016): “Advances in geometrical parametrization and reduced order models and methods for computational fluid dynamics problems in applied sciences and engineering: Overview and perspectives,” ECCOMAS Congress 2016 - Proceedings of the 7th European Congress on Computational Methods in Applied Sciences and Engineering, 1, 1013–1031. [Salmoiraghi et al.(2018)Salmoiraghi, Scardigli, Telib, and Rozza] Salmoiraghi, F., A. Scardigli, H. Telib, and G. Rozza (2018): “Free-form deformation, mesh morphing and reduced-order methods: enablers for efficient aerody- namic shape optimisation,” International Journal of Computational Fluid Dynamics, 32, 233–247. [Sandwell(1987)] Sandwell, D. T. (1987): “Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data,” Geophysical research letters, 14, 139–142. [Schilders et al.(2008)Schilders, van der Vorst, and Rommes] Schilders, W. H., H. A. van der Vorst, and J. Rommes (2008): Model Order Reduction: Theory, Research Aspects and Applications, Springer-Verlag Berlin Heidelberg. [Sederberg and Parry(1986)] Sederberg, T. and S. Parry (1986): “Free-Form Deformation of solid geometric models,” in Proceedings of SIGGRAPH
- Special Interest Group on GRAPHics and Interactive Techniques, SIG- GRAPH, 151–159. [Shepard(1968)] Shepard, D. (1968): “A two-dimensional interpolation function for irregularly-spaced data,” in Proceedings-1968 ACM National Conference, ACM, 517–524. [Sieger et al.(2015)Sieger, Menzel, and Botsch] Sieger, D., S. Menzel, and M. Botsch (2015): “On shape deformation techniques for simulation-based design optimization,” in New Challenges in Grid Generation and Adaptivity for Scientific Computing, Springer, 281–303.
Page 55
REFERENCES 55 [Smith(1985)] Smith, G. D. (1985): Numerical Solution of Partial Differential Equations: Finite Difference Methods, Clarendon Press, Oxford Applied Mathematics and Computing Science Series. [Tenenbaum et al.(2000)Tenenbaum, De Silva, and Langford] Tenenbaum, J. B., V. De Silva, and J. C. Langford (2000): “A global geometric framework for nonlinear dimensionality reduction,” science, 290, 2319–2323. [Tezzele et al.(2018a)Tezzele, Ballarin, and Rozza] Tezzele, M., F. Ballarin, and G. Rozza (2018a): “Combined parameter and model reduction of cardiovas- cular problems by means of active subspaces and POD-Galerkin methods,” in D. Boffi, L. F. Pavarino, G. Rozza, S. Scacchi, and C. Vergara, eds., Mathematical and Numerical Modeling of the Cardiovascular System and Applications, Springer International Publishing, 185–207. [Tezzele et al.(2018b)Tezzele, Demo, Gadalla, Mola, and Rozza] Tezzele, M., N. Demo, M. Gadalla, A. Mola, and G. Rozza (2018b): “Model order reduction by means of active subspaces and dynamic mode decomposition for parametric hull shape design hydrodynamics,” in Technology and Science for the Ships of the Future: Proceedings of NAV 2018: 19th International Conference on Ship & Maritime Research, IOS Press, 569–576. [Tezzele et al.(2019)Tezzele, Demo, and Rozza] Tezzele, M., N. Demo, and G. Rozza (2019): “Shape optimization through proper orthogonal decom- position with interpolation and dynamic mode decomposition enhanced by active subspaces,” in Proceedings of MARINE 2019: VIII International Conference on Computational Methods in Marine Engineering, 122–133. [Tezzele et al.(2018c)Tezzele, Salmoiraghi, Mola, and Rozza] Tezzele, M., F. Salmoiraghi, A. Mola, and G. Rozza (2018c): “Dimension reduction in heterogeneous parametric spaces with application to naval engineering shape design problems,” Advanced Modeling and Simulation in Engineering Sciences, 5, 25. [Vallaghé et al.(2011)Vallaghé, Fouquembergh, Le Hyaric, and Prud’Homme] Vallaghé, S., M. Fouquembergh, A. Le Hyaric, and C. Prud’Homme (2011): “A successive constraint method with minimal offline constraints for lower bounds of parametric coercivity constant,” https://hal.archives-ouvertes.fr/hal-00609212. [Van Der Maaten et al.(2009)Van Der Maaten, Postma, and Van den Herik] Van Der Maaten, L., E. Postma, and J. Van den Herik (2009): “Di- mensionality reduction: a comparative review,” J Mach Learn Res, 10, 66–71. [Verfuerth(2013)] Verfuerth, R. (2013): A posteriori error estimation techniques for finite element methods, Oxford Univ. Press.
Page 56
56 REFERENCES [Veroy and Patera(2005)] Veroy, K. and A. T. Patera (2005): “Certified real- time solution of the parametrized steady incompressible Navier–Stokes equations: rigorous reduced-basis a posteriori error bounds,” International Journal for Numerical Methods in Fluids, 47, 773–788. [Veroy et al.(2002)Veroy, Rovas, and Patera] Veroy, K., D. V. Rovas, and A. T. Patera (2002): “A posteriori error estimation for reduced-basis approxima- tion of parametrized elliptic coercive partial differential equations : “convex inverse” bound conditioners,” ESAIM: Control, Optimisation and Calculus of Variations, 8, 1007–1028. [Witteveen and Bijl(2009)] Witteveen, J. and H. Bijl (2009): “Explicit mesh deformation using Inverse Distance Weighting interpolation,” in 19th AIAA Computational Fluid Dynamics, AIAA. [Yano(2014)] Yano, M. (2014): “A space-time Petrov–Galerkin certified reduced basis method: Application to the boussinesq equations,” SIAM Journal on Scientific Computing, 36, A232–A266. [Yosida(1995)] Yosida, K. (1995): Functional Analysis, Springer-Verlag Berlin Heidelberg. [Zhang(2011)] Zhang, S. (2011): “Efficient greedy algorithms for successive con- straints methods with high-dimensional parameters,” Tech. Report 23, http://www.dam.brown.edu/people/shzhang/greedy_scm.pdf.
Canonical Hub: CANONICAL_INDEX